Ignore:
Timestamp:
Mar 30, 2011 9:31:40 AM (11 years ago)
Author:
raasch
Message:

formatting adjustments

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/advec_ws.f90

    r706 r709  
    11 MODULE advec_ws
    2 
    32
    43!-----------------------------------------------------------------------------!
    54! Current revisions:
    65! -----------------
     6! formatting adjustments
    77!
    88! Former revisions:
     
    4444!-----------------------------------------------------------------------------!
    4545
    46 
    4746    PRIVATE
    4847    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, &
     
    5251       MODULE PROCEDURE ws_init
    5352    END INTERFACE ws_init
     53
    5454    INTERFACE ws_statistics
    5555       MODULE PROCEDURE ws_statistics
    5656    END INTERFACE ws_statistics
     57
    5758    INTERFACE advec_s_ws
    5859       MODULE PROCEDURE advec_s_ws
    5960       MODULE PROCEDURE advec_s_ws_ij
    6061    END INTERFACE advec_s_ws
     62
    6163    INTERFACE advec_u_ws
    6264       MODULE PROCEDURE advec_u_ws
    6365       MODULE PROCEDURE advec_u_ws_ij
    6466    END INTERFACE advec_u_ws
     67
    6568    INTERFACE advec_v_ws
    6669       MODULE PROCEDURE advec_v_ws
    6770       MODULE PROCEDURE advec_v_ws_ij
    6871    END INTERFACE advec_v_ws
     72
    6973    INTERFACE advec_w_ws
    7074       MODULE PROCEDURE advec_w_ws
    7175       MODULE PROCEDURE advec_w_ws_ij
    7276    END INTERFACE advec_w_ws
     77
    7378    INTERFACE local_diss
    7479       MODULE PROCEDURE local_diss
     
    7883 CONTAINS
    7984
    80 !
    81 !-- Initialisation
    82 
     85
     86!------------------------------------------------------------------------------!
     87! Initialization of WS-scheme
     88!------------------------------------------------------------------------------!
    8389    SUBROUTINE ws_init
     90
    8491       USE arrays_3d
    8592       USE constants
     
    9198!--    Allocate arrays needed for dissipation control.
    9299       IF ( dissipation_control )  THEN
    93    !          ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr),  &
    94     !                  var_y(nzb+1:nzt,nys:nyn,nxl:nxr),  &
    95      !                 var_z(nzb+1:nzt,nys:nyn,nxl:nxr),  &
    96       !                gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    97        !               gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    98         !              gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
     100!          ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr),        &
     101!                   var_y(nzb+1:nzt,nys:nyn,nxl:nxr),        &
     102!                   var_z(nzb+1:nzt,nys:nyn,nxl:nxr),        &
     103!                   gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     104!                   gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     105!                   gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
    99106       ENDIF
    100              
     107
     108!
    101109!--    Set the appropriate factors for scalar and momentum advection.
    102110       adv_sca_5 = 0.016666666666666
    103        adv_sca_3 = 0.083333333333333 
     111       adv_sca_3 = 0.083333333333333
    104112       adv_mom_5 = 0.0083333333333333
    105113       adv_mom_3 = 0.041666666666666
     
    108116!--    Arrays needed for statical evaluation of fluxes.
    109117       IF ( ws_scheme_mom )  THEN
    110           ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:statistic_regions),         &
    111               sums_wsvs_ws_l(nzb:nzt+1,0:statistic_regions),               &
    112               sums_us2_ws_l(nzb:nzt+1,0:statistic_regions),                &
    113               sums_vs2_ws_l(nzb:nzt+1,0:statistic_regions),                &
    114               sums_ws2_ws_l(nzb:nzt+1,0:statistic_regions))
    115 
     118
     119          ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:statistic_regions),  &
     120                    sums_wsvs_ws_l(nzb:nzt+1,0:statistic_regions),  &
     121                    sums_us2_ws_l(nzb:nzt+1,0:statistic_regions),   &
     122                    sums_vs2_ws_l(nzb:nzt+1,0:statistic_regions),   &
     123                    sums_ws2_ws_l(nzb:nzt+1,0:statistic_regions))
     124
     125          sums_wsus_ws_l = 0.0
     126          sums_wsvs_ws_l = 0.0
     127          sums_us2_ws_l  = 0.0
     128          sums_vs2_ws_l  = 0.0
     129          sums_ws2_ws_l  = 0.0
     130
     131       ENDIF
     132
     133       IF ( ws_scheme_sca )  THEN
     134
     135          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:statistic_regions) )
     136          sums_wspts_ws_l = 0.0
     137
     138          IF ( humidity .OR. passive_scalar )  THEN
     139             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:statistic_regions) )
     140             sums_wsqs_ws_l = 0.0
     141          ENDIF
     142
     143          IF ( ocean )  THEN
     144             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:statistic_regions) )
     145             sums_wssas_ws_l = 0.0
     146          ENDIF
     147
     148       ENDIF
     149
     150!
     151!--    Arrays needed for reasons of speed optimization for cache and noopt
     152!--    version. For the vector version the buffer arrays are not necessary,
     153!--    because the the fluxes can swapped directly inside the loops of the
     154!--    advection routines.
     155       IF ( loop_optimization /= 'vector' )  THEN
     156
     157          IF ( ws_scheme_mom )  THEN
     158
     159             ALLOCATE( flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt),  &
     160                       flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt),  &
     161                       diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt) )
     162             ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn),               &
     163                       flux_l_v(nzb+1:nzt,nys:nyn),               &
     164                       flux_l_w(nzb+1:nzt,nys:nyn),               &
     165                       diss_l_u(nzb+1:nzt,nys:nyn),               &
     166                       diss_l_v(nzb+1:nzt,nys:nyn),               &
     167                       diss_l_w(nzb+1:nzt,nys:nyn) )
     168
     169          ENDIF
     170
     171          IF ( ws_scheme_sca )  THEN
     172
     173             ALLOCATE( flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt), &
     174                       diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt) )
     175             ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn),              &
     176                       flux_l_e(nzb+1:nzt,nys:nyn),               &
     177                       diss_l_pt(nzb+1:nzt,nys:nyn),              &
     178                       diss_l_e(nzb+1:nzt,nys:nyn) )
     179
     180             IF ( humidity .OR. passive_scalar )  THEN
     181                ALLOCATE( flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt) )
     182                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn),            &
     183                          diss_l_q(nzb+1:nzt,nys:nyn) )
     184             ENDIF
     185
     186             IF ( ocean )  THEN
     187                ALLOCATE( flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt) )
     188                ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn),           &
     189                          diss_l_sa(nzb+1:nzt,nys:nyn) )
     190             ENDIF
     191
     192          ENDIF
     193
     194       ENDIF
     195
     196!
     197!--    Determine the flags where the order of the scheme for horizontal
     198!--    advection has to be degraded.
     199       ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
     200       DO  i = nxl, nxr
     201          DO  j = nys, nyn
     202
     203             boundary_flags(j,i) = 0
     204             IF ( outflow_l )  THEN
     205                IF ( i == nxlu  )  boundary_flags(j,i) = 5
     206                IF ( i == nxl   )  boundary_flags(j,i) = 6
     207             ELSEIF ( outflow_r )  THEN
     208                IF ( i == nxr-1 )  boundary_flags(j,i) = 1
     209                IF ( i == nxr   )  boundary_flags(j,i) = 2
     210             ELSEIF ( outflow_n )  THEN
     211                IF ( j == nyn-1 )  boundary_flags(j,i) = 3
     212                IF ( j == nyn   )  boundary_flags(j,i) = 4
     213             ELSEIF ( outflow_s )  THEN
     214                IF ( j == nysv  )  boundary_flags(j,i) = 7
     215                IF ( j == nys   )  boundary_flags(j,i) = 8
     216             ENDIF
     217
     218          ENDDO
     219       ENDDO
     220       
     221    END SUBROUTINE ws_init
     222
     223
     224!------------------------------------------------------------------------------!
     225! Initialize variables used for storing statistic qauntities (fluxes, variances)
     226!------------------------------------------------------------------------------!
     227    SUBROUTINE ws_statistics
     228   
     229       USE control_parameters
     230       USE statistics
     231
     232       IMPLICIT NONE
     233
     234!       
     235!--    The arrays needed for statistical evaluation are set to to 0 at the
     236!--    begin of prognostic_equations.
     237       IF ( ws_scheme_mom )  THEN
    116238          sums_wsus_ws_l = 0.0
    117239          sums_wsvs_ws_l = 0.0
     
    122244
    123245       IF ( ws_scheme_sca )  THEN
    124           ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:statistic_regions))
    125246          sums_wspts_ws_l = 0.0
    126           IF ( humidity  .OR.  passive_scalar )  THEN
    127              ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:statistic_regions))
    128              sums_wsqs_ws_l = 0.0
    129           ENDIF
    130           IF ( ocean )  THEN
    131              ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:statistic_regions))
    132              sums_wssas_ws_l = 0.0
    133           ENDIF
     247          IF ( humidity .OR. passive_scalar )  sums_wsqs_ws_l = 0.0
     248          IF ( ocean )  sums_wssas_ws_l = 0.0
     249
    134250       ENDIF
    135 !
    136 !--    Arrays needed for reasons of speed optimization for cache and noopt
    137 !--    version. For the vector version the buffer arrays are not necessary,
    138 !--    because the the fluxes can swapped directly inside the loops of the
    139 !--    advection routines.
    140        IF ( loop_optimization /= 'vector' )  THEN
    141           IF ( ws_scheme_mom )  THEN
    142              ALLOCATE(flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt),            &
    143                 flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt),                  &
    144                 diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt))
    145              ALLOCATE(flux_l_u(nzb+1:nzt,nys:nyn),                         &
    146                 flux_l_v(nzb+1:nzt,nys:nyn), flux_l_w(nzb+1:nzt,nys:nyn),  &
    147                 diss_l_u(nzb+1:nzt,nys:nyn), diss_l_v(nzb+1:nzt,nys:nyn),  &
    148                 diss_l_w(nzb+1:nzt,nys:nyn))
    149           ENDIF
    150           IF ( ws_scheme_sca )  THEN
    151              ALLOCATE(flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt),           &
    152                 diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt))
    153              ALLOCATE(flux_l_pt(nzb+1:nzt,nys:nyn),                        &
    154                 flux_l_e(nzb+1:nzt,nys:nyn), diss_l_pt(nzb+1:nzt,nys:nyn), &
    155                 diss_l_e(nzb+1:nzt,nys:nyn))
    156              IF ( humidity  .OR.  passive_scalar )  THEN
    157                 ALLOCATE(flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt))
    158                 ALLOCATE(flux_l_q(nzb+1:nzt,nys:nyn),                      &
    159                    diss_l_q(nzb+1:nzt,nys:nyn))
    160              ENDIF
    161              IF ( ocean )  THEN
    162                 ALLOCATE(flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt))
    163                 ALLOCATE(flux_l_sa(nzb+1:nzt,nys:nyn),                     &
    164                    diss_l_sa(nzb+1:nzt,nys:nyn))
    165              ENDIF
    166           ENDIF
    167        ENDIF
    168 !
    169 !--    Determine the flags where the order of the scheme for horizontal
    170 !--    advection should be degraded.
    171        ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
    172        DO  i = nxl, nxr
    173           DO  j = nys, nyn
    174              boundary_flags(j,i) = 0
    175              IF ( outflow_l )  THEN
    176                 IF ( i == nxlu )  boundary_flags(j,i) = 5
    177                 IF ( i == nxl )  boundary_flags(j,i) = 6
    178              ELSEIF ( outflow_r )  THEN
    179                 IF ( i == nxr-1 )  boundary_flags(j,i) = 1
    180                 IF ( i == nxr )  boundary_flags(j,i) = 2
    181              ELSEIF ( outflow_n )  THEN
    182                 IF ( j == nyn-1 )  boundary_flags(j,i) = 3
    183                 IF ( j == nyn )  boundary_flags(j,i) = 4
    184              ELSEIF ( outflow_s )  THEN
    185                 IF ( j == nysv )  boundary_flags(j,i) = 7
    186                 IF ( j == nys )  boundary_flags(j,i) = 8
    187              ENDIF
    188           ENDDO
    189        ENDDO
    190        
    191     END SUBROUTINE ws_init
    192    
    193 
    194    
    195     SUBROUTINE ws_statistics
    196    
    197        USE control_parameters
    198        USE statistics
    199 
    200 !       
    201 !--      The arrays needed for statistical evaluation are set to to 0 at the
    202 !--      begin of prognostic_equations.
    203           IF ( ws_scheme_mom )  THEN
    204              sums_wsus_ws_l = 0.0
    205              sums_wsvs_ws_l = 0.0
    206              sums_us2_ws_l = 0.0
    207              sums_vs2_ws_l = 0.0
    208              sums_ws2_ws_l = 0.0
    209           ENDIF
    210           IF ( ws_scheme_sca )  THEN
    211              sums_wspts_ws_l = 0.0
    212              IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0
    213              IF ( ocean )  sums_wssas_ws_l = 0.0
    214           ENDIF
    215251
    216252    END SUBROUTINE ws_statistics
    217 
    218253
    219254
     
    222257!------------------------------------------------------------------------------!
    223258    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
    224                swap_diss_y_local, swap_flux_x_local, swap_diss_x_local  )
     259                              swap_diss_y_local, swap_flux_x_local, &
     260                              swap_diss_x_local  )
    225261
    226262       USE arrays_3d
     
    239275       REAL, DIMENSION(nzb:nzt+1)         :: flux_t, diss_t, flux_r, diss_r,  &
    240276                                             flux_n, diss_n
    241        REAL, DIMENSION(nzb+1:nzt)         :: swap_flux_y_local,               &       
     277       REAL, DIMENSION(nzb+1:nzt)         :: swap_flux_y_local,               &
    242278                                             swap_diss_y_local
    243279       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local,               &
    244280                                             swap_diss_x_local
    245281       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
    246        
     282
     283
    247284       degraded_l = .FALSE.
    248285       degraded_s = .FALSE.
     
    250287       IF ( boundary_flags(j,i) /= 0 )  THEN
    251288!       
    252 !--      Degrade the order for Dirichlet bc. at the outflow boundary.
    253          SELECT CASE ( boundary_flags(j,i) )
    254          
    255             CASE ( 1 )
    256                DO  k = nzb_s_inner(j,i) + 1, nzt
    257                   u_comp    = u(k,j,i+1) - u_gtrans
    258                   flux_r(k) = u_comp * (                                      &
    259                               7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
    260                               -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    261                   diss_r(k) = - abs(u_comp) * (                               &
    262                               3. * ( sk(k,j,i+1)-sk(k,j,i)     )              &
    263                               -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     289!--       Degrade the order for Dirichlet bc. at the outflow boundary
     290          SELECT CASE ( boundary_flags(j,i) )
     291
     292             CASE ( 1 )
     293
     294                DO  k = nzb_s_inner(j,i)+1, nzt
     295                   u_comp    = u(k,j,i+1) - u_gtrans
     296                   flux_r(k) = u_comp * (                                      &
     297                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
     298                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
     299                   diss_r(k) = -ABS( u_comp ) * (                              &
     300                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
     301                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     302                ENDDO
     303
     304             CASE ( 2 )
     305
     306                DO  k = nzb_s_inner(j,i)+1, nzt
     307                   u_comp    = u(k,j,i+1) - u_gtrans
     308                   flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
     309                   diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i),  &
     310                                         sk(k,j,i-1), sk(k,j,i-2), u_comp,     &
     311                                         0.5, ddx )
     312                ENDDO
     313
     314             CASE ( 3 )
     315
     316                DO  k = nzb_s_inner(j,i) + 1, nzt
     317                   v_comp = v(k,j+1,i) - v_gtrans
     318                   flux_n(k) = v_comp * (                                      &
     319                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
     320                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
     321                   diss_n(k) = -ABS( v_comp ) * (                              &
     322                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)     )           &
     323                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    264324               ENDDO
    265                
    266             CASE ( 2 )
    267                DO  k = nzb_s_inner(j,i) + 1, nzt
    268                   u_comp    = u(k,j,i+1) - u_gtrans
    269                   flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
    270                   diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i),  &
    271                                         sk(k,j,i-1), sk(k,j,i-2), u_comp,     &
    272                                         0.5, ddx )
    273               ENDDO
    274              
    275             CASE ( 3 )
    276                DO  k = nzb_s_inner(j,i) + 1, nzt
    277                   v_comp = v(k,j+1,i) - v_gtrans
    278                   flux_n(k) = v_comp * (                                      &
    279                               7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
    280                               -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    281                   diss_n(k) = - abs(v_comp) * (                               &
    282                               3. * ( sk(k,j+1,i)-sk(k,j,i)    )               &
    283                               -    (sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    284               ENDDO
    285             CASE ( 4 )
    286                DO  k = nzb_s_inner(j,i) + 1, nzt
    287                   v_comp    = v(k,j+1,i) - v_gtrans
    288                   flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
    289                   diss_n(k) =  diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
     325
     326             CASE ( 4 )
     327
     328                DO  k = nzb_s_inner(j,i)+1, nzt
     329                   v_comp    = v(k,j+1,i) - v_gtrans
     330                   flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
     331                   diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
    290332                                         sk(k,j-1,i), sk(k,j-2,i), v_comp,    &
    291333                                         0.5, ddy )     
    292                ENDDO
     334                ENDDO
     335
     336             CASE ( 5 )
     337
     338                DO  k = nzb_w_inner(j,i)+1, nzt
     339                   u_comp    = u(k,j,i+1) - u_gtrans
     340                   flux_r(k) = u_comp  * (                                     &
     341                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
     342                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
     343                   diss_r(k) = -ABS( u_comp ) * (                              &
     344                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
     345                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     346                ENDDO 
     347
     348             CASE ( 6 )
     349
     350                DO  k = nzb_s_inner(j,i)+1, nzt
     351                   u_comp    = u(k,j,i+1) - u_gtrans
     352                   flux_r(k) = u_comp * (                                      &
     353                               7.0 * ( sk(k,j,i+1) + sk(k,j,i)   )             &
     354                                   - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
     355                   diss_r(k) = -ABS( u_comp ) * (                              &
     356                               3.0 * ( sk(k,j,i+1) - sk(k,j,i)   )             &
     357                                   - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
     358!                             
     359!--                Compute leftside fluxes for the left boundary of PE domain
     360                   u_comp                 = u(k,j,i) - u_gtrans
     361                   swap_flux_x_local(k,j) = u_comp * (                         &
     362                                            sk(k,j,i) + sk(k,j,i-1) ) * 0.5
     363                   swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
     364                                                      sk(k,j,i), sk(k,j,i-1),  &
     365                                                      sk(k,j,i-1), u_comp,     &
     366                                                      0.5, ddx )
     367                ENDDO
     368                degraded_l = .TRUE.
     369
     370             CASE ( 7 )
     371
     372                DO  k = nzb_s_inner(j,i)+1, nzt
     373                   v_comp    = v(k,j+1,i)-v_gtrans
     374                   flux_n(k) = v_comp * (                                      &
     375                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
     376                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
     377                   diss_n(k) = -ABS( v_comp ) * (                              &
     378                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
     379                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
     380                ENDDO
     381
     382             CASE ( 8 )
     383
     384                DO  k = nzb_s_inner(j,i)+1, nzt
     385                   v_comp    = v(k,j+1,i) - v_gtrans
     386                   flux_n(k) = v_comp * (                                      &
     387                               7.0 * ( sk(k,j+1,i) + sk(k,j,i)   )             &
     388                                   - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
     389                   diss_n(k) = -ABS( v_comp ) * (                              &
     390                               3.0 * ( sk(k,j+1,i) - sk(k,j,i)   )             &
     391                                   - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
     392!                             
     393!--                Compute southside fluxes for the south boundary of PE domain
     394                   v_comp               = v(k,j,i) - v_gtrans
     395                   swap_flux_y_local(k) = v_comp *                             &
     396                                          ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5
     397                   swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
     398                                                    sk(k,j,i), sk(k,j-1,i),    &
     399                                                    sk(k,j-1,i), v_comp,       &
     400                                                    0.5, ddy )
     401                ENDDO
     402                degraded_s = .TRUE.
    293403               
    294             CASE ( 5 )
    295                DO  k = nzb_w_inner(j,i) + 1, nzt
    296                   u_comp    = u(k,j,i+1) - u_gtrans
    297                   flux_r(k) = u_comp  * (                                     &
    298                               7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
    299                               -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    300                   diss_r(k) = - abs(u_comp) * (                               &
    301                               3. * ( sk(k,j,i+1)-sk(k,j,i) )                  &
    302                               -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
    303                ENDDO 
    304                
    305             CASE ( 6 )
    306                DO  k = nzb_s_inner(j,i) + 1, nzt
    307                   u_comp    = u(k,j,i+1) - u_gtrans
    308                   flux_r(k) = u_comp * (                                      &
    309                               7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
    310                               -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
    311                   diss_r(k) = - abs(u_comp) * (                               &
    312                               3. * ( sk(k,j,i+1)-sk(k,j,i)     )              &
    313                               -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
    314 !                             
    315 !--               Compute leftside fluxes for the left boundary of PE domain.
    316                   u_comp                 = u(k,j,i) - u_gtrans
    317                   swap_flux_x_local(k,j) = u_comp * (                         &
    318                                            sk(k,j,i) + sk(k,j,i-1) ) * 0.5
    319                   swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
    320                                                      sk(k,j,i), sk(k,j,i-1),  &
    321                                                      sk(k,j,i-1), u_comp,     &
    322                                                      0.5, ddx )
    323                ENDDO
    324                degraded_l = .TRUE.
    325                
    326             CASE ( 7 )
    327                DO  k = nzb_s_inner(j,i) + 1, nzt
    328                   v_comp    = v(k,j+1,i)-v_gtrans
    329                   flux_n(k) = v_comp * (                                      &
    330                               7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
    331                               -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    332                   diss_n(k) = - abs(v_comp) * (                               &
    333                               3. * ( sk(k,j+1,i) - sk(k,j,i)   )              &
    334                               -    ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    335                ENDDO
    336                
    337             CASE ( 8 )
    338                DO  k = nzb_s_inner(j,i) + 1, nzt
    339                   v_comp    = v(k,j+1,i) - v_gtrans
    340                   flux_n(k) = v_comp * (                                      &
    341                               7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
    342                               -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
    343                   diss_n(k) = - abs(v_comp) * (                               &
    344                               3. * ( sk(k,j+1,i) - sk(k,j,i) )                &
    345                               -    ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
    346 !                             
    347 !--               Compute southside fluxes for the south boundary of PE domain           
    348                   v_comp               = v(k,j,i) - v_gtrans
    349                   swap_flux_y_local(k) = v_comp * (                           &
    350                                          sk(k,j,i)+ sk(k,j-1,i) ) * 0.5
    351                   swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
    352                                                    sk(k,j,i), sk(k,j-1,i),    &
    353                                                    sk(k,j-1,i), v_comp,       &
    354                                                    0.5, ddy )
    355                ENDDO
    356                degraded_s = .TRUE.
    357                
    358             CASE DEFAULT
     404             CASE DEFAULT
    359405           
    360          END SELECT
     406          END SELECT
     407
    361408!         
    362 !--      Compute the crosswise 5th order fluxes at the outflow
    363          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
    364          .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
    365             DO  k = nzb_s_inner(j,i) + 1, nzt
    366                v_comp    = v(k,j+1,i) - v_gtrans
    367                flux_n(k) = v_comp * (                                         &
    368                            37.  * ( sk(k,j+1,i) + sk(k,j,i)   )               &
    369                            - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )               &
    370                            +      ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
    371                diss_n(k) = - abs(v_comp) * (                                  &
    372                            10. *  ( sk(k,j+1,i) - sk(k,j,i)   )               &
    373                            - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )               &
    374                            +      ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
    375             ENDDO
     409!--       Compute the crosswise 5th order fluxes at the outflow
     410          IF ( boundary_flags(j,i) == 1  .OR.  boundary_flags(j,i) == 2  .OR. &
     411               boundary_flags(j,i) == 5  .OR.  boundary_flags(j,i) == 6 )  THEN
     412
     413             DO  k = nzb_s_inner(j,i)+1, nzt
     414                v_comp    = v(k,j+1,i) - v_gtrans
     415                flux_n(k) = v_comp * (                                         &
     416                            37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )               &
     417                          -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )               &
     418                          +        ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
     419                diss_n(k) = -ABS( v_comp ) * (                                 &
     420                            10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )               &
     421                          -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )               &
     422                          +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
     423             ENDDO
    376424           
    377          ELSE
     425          ELSE
    378426         
    379             DO  k = nzb_s_inner(j,i) + 1, nzt
    380                u_comp    = u(k,j,i+1) - u_gtrans 
    381                flux_r(k) = u_comp * (                                         &
    382                            37. * ( sk(k,j,i+1) + sk(k,j,i)   )               &
    383                            - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )               &
    384                            +      ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
    385                diss_r(k) = - abs(u_comp) * (                                  &
    386                            10. * ( sk(k,j,i+1) - sk(k,j,i)   )               &
    387                            - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )               &
    388                            +      ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    389             ENDDO
     427             DO  k = nzb_s_inner(j,i)+1, nzt
     428                u_comp    = u(k,j,i+1) - u_gtrans 
     429                flux_r(k) = u_comp * (                                         &
     430                            37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )               &
     431                          -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )               &
     432                                ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
     433                diss_r(k) = -ABS( u_comp ) * (                                 &
     434                            10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )               &
     435                          -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )               &
     436                                ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
     437             ENDDO
    390438           
    391          ENDIF
    392 !
    393 !--    Compute the fifth order fluxes for the interior of PE domain.
     439          ENDIF
     440
    394441       ELSE
    395442               
    396           DO  k = nzb_u_inner(j,i) + 1, nzt
     443!
     444!--       Compute the fifth order fluxes for the interior of PE domain.
     445          DO  k = nzb_u_inner(j,i)+1, nzt
    397446             u_comp    = u(k,j,i+1) - u_gtrans
    398447             flux_r(k) = u_comp * (                                           &
    399                          37.  * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
    400                          - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
    401                          +      ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
    402              diss_r(k) = - abs(u_comp) * (                                    &
    403                          10.  * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
    404                          - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
    405                          +      ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
     448                         37.0 * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
     449                       -  8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
     450                             ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
     451             diss_r(k) = -ABS( u_comp ) * (                                   &
     452                         10.0 * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
     453                       -  5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
     454                             ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
    406455
    407456             v_comp    = v(k,j+1,i) - v_gtrans
    408457             flux_n(k) = v_comp * (                                           &
    409                          37.  * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
    410                          - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
    411                          +      ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
    412              diss_n(k) = - abs(v_comp) * (                                    &
    413                          10.  * ( sk(k,j+1,i) - sk(k,j,i)   )                 &                     
    414                          - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
    415                          +      ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5               
     458                         37.0 * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
     459                       -  8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
     460                             ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
     461             diss_n(k) = -ABS( v_comp ) * (                                   &
     462                         10.0 * ( sk(k,j+1,i) - sk(k,j,i)   )                 &
     463                       -  5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
     464                       +        ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
    416465          ENDDO
    417466         
     
    424473              v_comp               = v(k,j,i) - v_gtrans
    425474              swap_flux_y_local(k) = v_comp * (                               &
    426                                      37. *    ( sk(k,j,i) + sk(k,j-1,i)   )   &
    427                                      - 8. *   ( sk(k,j+1,i) + sk(k,j-2,i) )   &
    428                                      +        ( sk(k,j+2,i) + sk(k,j-3,i) )   &
    429                                      ) * adv_sca_5
    430               swap_diss_y_local(k) = - abs(v_comp) * (                        &
    431                                      10. *  ( sk(k,j,i) - sk(k,j-1,i)  )     &
    432                                      - 5. * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
    433                                      +        sk(k,j+2,i) - sk(k,j-3,i)       &
    434                                      ) * adv_sca_5
     475                                     37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
     476                                   -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
     477                                   +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
     478                                              ) * adv_sca_5
     479              swap_diss_y_local(k) = -ABS( v_comp ) * (                       &
     480                                     10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
     481                                   -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
     482                                           sk(k,j+2,i) - sk(k,j-3,i)       &
     483                                                      ) * adv_sca_5
    435484           ENDDO           
    436485         
    437486        ENDIF
    438487       
    439         IF ( i == nxl .AND. ( .NOT. degraded_l ) )  THEN
     488        IF ( i == nxl  .AND.  .NOT. degraded_l )  THEN
    440489       
    441            DO  k = nzb_s_inner(j,i) + 1, nzt
     490           DO  k = nzb_s_inner(j,i)+1, nzt
    442491              u_comp                 = u(k,j,i) - u_gtrans
    443492              swap_flux_x_local(k,j) = u_comp * (                             &
    444                                        37.  * ( sk(k,j,i) + sk(k,j,i-1)  )   &
    445                                        - 8. * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
    446                                        +      ( sk(k,j,i+2) + sk(k,j,i-3) )   &
    447                                        ) * adv_sca_5
    448               swap_diss_x_local(k,j) = - abs(u_comp) * (                      &
    449                                        10.  * ( sk(k,j,i) - sk(k,j,i-1)  )   &
    450                                        - 5. * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
    451                                        +      ( sk(k,j,i+2) - sk(k,j,i-3) )   &
    452                                        ) * adv_sca_5
     493                                       37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
     494                                     -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
     495                                           ( sk(k,j,i+2) + sk(k,j,i-3) )   &
     496                                                ) * adv_sca_5
     497              swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
     498                                       10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
     499                                     -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
     500                                           ( sk(k,j,i+2) - sk(k,j,i-3) )   &
     501                                                        ) * adv_sca_5
    453502           ENDDO
    454503           
    455504        ENDIF
     505
    456506!       
    457 !--    Now compute the tendency terms for the horizontal parts.
    458        DO  k = nzb_s_inner(j,i) + 1, nzt
     507!--    Now compute the tendency terms for the horizontal parts
     508       DO  k = nzb_s_inner(j,i)+1, nzt
    459509         
    460           tend(k,j,i) = tend(k,j,i) - (                                       &
    461                ( flux_r(k) + diss_r(k)                                        &
    462              - ( swap_flux_x_local(k,j) + swap_diss_x_local(k,j) ) ) * ddx    &
    463              + (flux_n(k)+diss_n(k)                                           &
    464              - (swap_flux_y_local(k)    + swap_diss_y_local(k) ) )   * ddy )
     510          tend(k,j,i) = tend(k,j,i) - (                                      &
     511                          ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
     512                            swap_diss_x_local(k,j) ) * ddx                   &
     513                        + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   - &
     514                            swap_diss_y_local(k)   ) * ddy                   &
     515                                      )
    465516
    466517          swap_flux_y_local(k)   = flux_n(k)
     
    470521         
    471522       ENDDO
    472 !
    473 !--    Vertical advection, degradation of order near surface and top.
    474 !--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
    475 !--    statistical evaluation the top flux at the surface should be 0
    476 
    477        flux_t(nzb_s_inner(j,i)) = 0.
    478        diss_t(nzb_s_inner(j,i)) = 0.
     523
     524!
     525!--    Vertical advection, degradation of order near bottom and top.
     526!--    The fluxes flux_d and diss_d at the surface are 0. Due to later
     527!--    calculation of statistics the top flux at the surface should be 0.
     528       flux_t(nzb_s_inner(j,i)) = 0.0
     529       diss_t(nzb_s_inner(j,i)) = 0.0
     530
    479531!       
    480 !--    2nd order scheme       
     532!--    2nd-order scheme (bottom)
    481533       k         = nzb_s_inner(j,i) + 1
    482534       flux_d    = flux_t(k-1)
    483535       diss_d    = diss_t(k-1)
    484536       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
     537
    485538!       
    486 !--    sk(k,j,i) is referenced three times to avoid a access below surface.
     539!--    sk(k,j,i) is referenced three times to avoid an access below surface
    487540       diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), sk(k,j,i),  &
    488541                             sk(k,j,i), w(k,j,i), 0.5, ddzw(k) )
    489542                   
    490        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                   &
    491                                  - ( flux_d + diss_d )     ) * ddzw(k)
    492 !
    493 !--    WS3 as an intermediate step      
     543       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
     544                                   * ddzw(k)
     545!
     546!--    WS3 as an intermediate step (bottom)
    494547       k         = nzb_s_inner(j,i) + 2
    495548       flux_d    = flux_t(k-1)
    496549       diss_d    = diss_t(k-1)
    497550       flux_t(k) = w(k,j,i) * (                                               &
    498                    7. * ( sk(k+1,j,i) + sk(k,j,i)   )                         &
    499                    -    ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
    500        diss_t(k) = - abs(w(k,j,i)) * (                                        &
    501                    3. * ( sk(k+1,j,i) - sk(k,j,i)   )                         &
    502                    -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
    503 
    504        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    505                                  - ( flux_d + diss_d ) ) * ddzw(k)
     551                                7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )           &
     552                                    - ( sk(k+2,j,i) + sk(k-1,j,i) )           &
     553                              ) * adv_sca_3
     554       diss_t(k) = -ABS( w(k,j,i) ) * (                                       &
     555                                3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )           &
     556                                    - ( sk(k+2,j,i) - sk(k-1,j,i) )           &
     557                                      ) * adv_sca_3
     558
     559       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
     560                                   * ddzw(k)
    506561!
    507562!--    WS5
    508        DO  k = nzb_s_inner(j,i) + 3, nzt - 2
     563       DO  k = nzb_s_inner(j,i)+3, nzt-2
    509564       
    510565          flux_d    = flux_t(k-1)
    511           diss_d    = diss_t(k-1)
    512          
    513           flux_t(k) = w(k,j,i) * (                                            & 
    514                       37.  * ( sk(k+1,j,i) + sk(k,j,i)   )                    &
    515                       - 8. * ( sk(k+2,j,i) + sk(k-1,j,i) )                    &
    516                       +      ( sk(k+3,j,i) + sk(k-2,j,i) ) ) *adv_sca_5
    517           diss_t(k) = - abs(w(k,j,i)) * (                                     &
    518                       10.  * ( sk(k+1,j,i) - sk(k,j,i)   )                    &     
    519                       - 5. * ( sk(k+2,j,i) - sk(k-1,j,i) )                    &
    520                       +      ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
     566          diss_d    = diss_t(k-1)         
     567          flux_t(k) = w(k,j,i) * (                                            &
     568                                   37.0 * ( sk(k+1,j,i) + sk(k,j,i)   )       &
     569                                 -  8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) )       &
     570                                 +        ( sk(k+3,j,i) + sk(k-2,j,i) )       &
     571                                 ) * adv_sca_5
     572          diss_t(k) = -ABS( w(k,j,i) ) * (                                     &
     573                                           10.0 * ( sk(k+1,j,i) - sk(k,j,i)   )&
     574                                         -  5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )&
     575                                         +        ( sk(k+3,j,i) - sk(k-2,j,i) )&
     576                                         ) * adv_sca_5
    521577
    522578          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
    523579                                    - ( flux_d + diss_d ) ) * ddzw(k)
     580
    524581       ENDDO
    525 !
    526 !--    WS3 as an intermediate step
     582
     583!
     584!--    WS3 as an intermediate step (top)
    527585       k         = nzt - 1
    528586       flux_d    = flux_t(k-1)
    529587       diss_d    = diss_t(k-1)
    530        flux_t(k) = w(k,j,i) * (                                               &
    531                    7. * ( sk(k+1,j,i) + sk(k,j,i) )                           &
    532                    -    ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
    533        diss_t(k) = - abs(w(k,j,i)) * (                                        &
    534                    3. * ( sk(k+1,j,i) - sk(k,j,i) )                           &
    535                    -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
    536 
    537        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
    538                                  - ( flux_d + diss_d ) ) * ddzw(k)
     588       flux_t(k) = w(k,j,i) * (                                                &
     589                                7.0 * ( sk(k+1,j,i) + sk(k,j,i)   )            &
     590                                    - ( sk(k+2,j,i) + sk(k-1,j,i) )            &
     591                              ) * adv_sca_3
     592       diss_t(k) = -ABS( w(k,j,i) ) * (                                        &
     593                                        3.0 * ( sk(k+1,j,i) - sk(k,j,i)   )    &
     594                                            - ( sk(k+2,j,i) - sk(k-1,j,i) )    &
     595                                      ) * adv_sca_3
     596
     597       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
     598                                   * ddzw(k)
    539599!                                 
    540 !--    2nd
     600!--    2nd-order scheme (top)
    541601       k         = nzt
    542602       flux_d    = flux_t(k-1)
    543603       diss_d    = diss_t(k-1)
    544        flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) *0.5
     604       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
     605
    545606!       
    546 !--    sk(k+1) is referenced two times to avoid a segmentation fault at top.
    547        diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i),&
     607!--    sk(k+1) is referenced two times to avoid a segmentation fault at top
     608       diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i), &
    548609                             sk(k-2,j,i), w(k,j,i), 0.5, ddzw(k) )
    549610
    550        tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                   &
    551                                  - ( flux_d + diss_d ) ) * ddzw(k)
     611       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
     612                                   * ddzw(k)
    552613!       
    553 !--    evaluation of statistics
     614!--    Evaluation of statistics
    554615       SELECT CASE ( sk_char )
    555616
    556617          CASE ( 'pt' )
    557             DO  k = nzb_s_inner(j,i), nzt
    558                sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:)                     &
    559                  + (flux_t(k)+diss_t(k))                                       &
    560                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     618
     619             DO  k = nzb_s_inner(j,i), nzt
     620                sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:) +                  &
     621                                       ( flux_t(k) + diss_t(k) )               &
     622                                 * weight_substep(intermediate_timestep_count) &
     623                                       * rmask(j,i,:)
    561624             ENDDO
    562625             
    563626          CASE ( 'sa' )
    564             DO  k = nzb_s_inner(j,i), nzt
    565                sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:)                     &
    566                  +(flux_t(k)+diss_t(k))                                        &
    567                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     627
     628             DO  k = nzb_s_inner(j,i), nzt
     629                sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:) +                  &
     630                                       ( flux_t(k) + diss_t(k) )               &
     631                                 * weight_substep(intermediate_timestep_count) &
     632                                       * rmask(j,i,:)
    568633             ENDDO
    569634             
    570635          CASE ( 'q' )
     636
    571637             DO  k = nzb_s_inner(j,i), nzt
    572                 sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:)                      &
    573                   + (flux_t(k)+diss_t(k))                                      &
    574                   * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     638                sums_wsqs_ws_l(k,:)  = sums_wsqs_ws_l(k,:) +                   &
     639                                      ( flux_t(k) + diss_t(k) )                &
     640                                 * weight_substep(intermediate_timestep_count) &
     641                                      * rmask(j,i,:)
    575642             ENDDO
    576643
Note: See TracChangeset for help on using the changeset viewer.