Ignore:
Timestamp:
Aug 17, 2011 2:13:26 PM (10 years ago)
Author:
suehring
Message:

Bugfix in ws-scheme concerning OpenMP paralellization.

File:
1 edited

Legend:

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

    r714 r736  
    33!-----------------------------------------------------------------------------!
    44! Current revisions:
    5 ! -----------------
     5! ------------------
     6! Bugfix concerning OpenMP parallelization. i_omp introduced, because first
     7! index where fluxes on left side have to be calculated explicitly is
     8! different on each thread. Furthermore the swapping of fluxes is now
     9! thread-safe by adding an additional dimension.
    610!
    711! Former revisions:
     
    4347! A further routine local_diss_ij is available (next weeks) and is used if a
    4448! control of dissipative fluxes is desired.
    45 ! In case of vertical grid stretching the order of advection scheme is
    46 ! degraded. This is also realized for the ocean where the stretching is
    47 ! downwards.
    4849!
    4950! OUTSTANDING: - Dirichlet and Radiation boundary conditions for
     
    102103       USE control_parameters
    103104       USE indices
     105       USE pegrid
    104106       USE statistics
    105107
     
    166168          IF ( ws_scheme_mom )  THEN
    167169
    168              ALLOCATE( flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt),  &
    169                        flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt),  &
    170                        diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt) )
    171              ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn),               &
    172                        flux_l_v(nzb+1:nzt,nys:nyn),               &
    173                        flux_l_w(nzb+1:nzt,nys:nyn),               &
    174                        diss_l_u(nzb+1:nzt,nys:nyn),               &
    175                        diss_l_v(nzb+1:nzt,nys:nyn),               &
    176                        diss_l_w(nzb+1:nzt,nys:nyn) )
     170             ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1),            &
     171                       flux_s_v(nzb+1:nzt,0:threads_per_task-1),            &
     172                       flux_s_w(nzb+1:nzt,0:threads_per_task-1),            &
     173                       diss_s_u(nzb+1:nzt,0:threads_per_task-1),            &
     174                       diss_s_v(nzb+1:nzt,0:threads_per_task-1),            &
     175                       diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
     176             ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
     177                       flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
     178                       flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
     179                       diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
     180                       diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
     181                       diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    177182
    178183          ENDIF
     
    180185          IF ( ws_scheme_sca )  THEN
    181186
    182              ALLOCATE( flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt), &
    183                        diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt) )
    184              ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn),              &
    185                        flux_l_e(nzb+1:nzt,nys:nyn),               &
    186                        diss_l_pt(nzb+1:nzt,nys:nyn),              &
    187                        diss_l_e(nzb+1:nzt,nys:nyn) )
     187             ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1),           &
     188                       flux_s_e(nzb+1:nzt,0:threads_per_task-1),            &
     189                       diss_s_pt(nzb+1:nzt,0:threads_per_task-1),           &
     190                       diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
     191             ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     192                       flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
     193                       diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     194                       diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    188195
    189196             IF ( humidity .OR. passive_scalar )  THEN
    190                 ALLOCATE( flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt) )
    191                 ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn),            &
    192                           diss_l_q(nzb+1:nzt,nys:nyn) )
     197                ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),         &
     198                          diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
     199                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
     200                          diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    193201             ENDIF
    194202
    195203             IF ( ocean )  THEN
    196                 ALLOCATE( flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt) )
    197                 ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn),           &
    198                           diss_l_sa(nzb+1:nzt,nys:nyn) )
     204                ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1),         &
     205                          diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
     206                ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
     207                          diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    199208             ENDIF
    200209
     
    267276    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
    268277                              swap_diss_y_local, swap_flux_x_local, &
    269                               swap_diss_x_local  )
     278                              swap_diss_x_local, i_omp, tn )
    270279
    271280       USE arrays_3d
     
    274283       USE grid_variables
    275284       USE indices
     285       USE pegrid
    276286       USE statistics
    277287
    278288       IMPLICIT NONE
    279289
    280        INTEGER ::  i, j, k
     290       INTEGER ::  i, i_omp, j, k, tn
    281291       LOGICAL ::  degraded_l, degraded_s
    282292       REAL    ::  flux_d, diss_d, u_comp, v_comp
     
    284294       REAL, DIMENSION(nzb:nzt+1)         :: flux_t, diss_t, flux_r, diss_r,  &
    285295                                             flux_n, diss_n
    286        REAL, DIMENSION(nzb+1:nzt)         :: swap_flux_y_local,               &
    287                                              swap_diss_y_local
    288        REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local,               &
    289                                              swap_diss_x_local
     296       REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local,  &
     297                                                          swap_diss_y_local
     298       REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ::             &
     299                                                          swap_flux_x_local,  &
     300                                                          swap_diss_x_local
    290301       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
    291302
     
    368379!--                Compute leftside fluxes for the left boundary of PE domain
    369380                   u_comp                 = u(k,j,i) - u_gtrans
    370                    swap_flux_x_local(k,j) = u_comp * (                         &
     381                   swap_flux_x_local(k,j,tn) = u_comp * (                         &
    371382                                            sk(k,j,i) + sk(k,j,i-1) ) * 0.5
    372                    swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
     383                   swap_diss_x_local(k,j,tn) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
    373384                                                      sk(k,j,i), sk(k,j,i-1),  &
    374385                                                      sk(k,j,i-1), u_comp,     &
     
    402413!--                Compute southside fluxes for the south boundary of PE domain
    403414                   v_comp               = v(k,j,i) - v_gtrans
    404                    swap_flux_y_local(k) = v_comp *                             &
     415                   swap_flux_y_local(k,tn) = v_comp *                             &
    405416                                          ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5
    406                    swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
     417                   swap_diss_y_local(k,tn) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
    407418                                                    sk(k,j,i), sk(k,j-1,i),    &
    408419                                                    sk(k,j-1,i), v_comp,       &
     
    477488!       
    478489!--    Compute left- and southside fluxes of the respective PE bounds.       
    479        IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
     490       IF ( j == nys  .AND.  .NOT. degraded_s )  THEN
    480491       
    481492           DO  k = nzb_s_inner(j,i)+1, nzt
    482493              v_comp               = v(k,j,i) - v_gtrans
    483               swap_flux_y_local(k) = v_comp * (                               &
     494              swap_flux_y_local(k,tn) = v_comp * (                               &
    484495                                     37.0 * ( sk(k,j,i)   + sk(k,j-1,i) )     &
    485496                                   -  8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) )     &
    486497                                   +        ( sk(k,j+2,i) + sk(k,j-3,i) )     &
    487498                                              ) * adv_sca_5
    488               swap_diss_y_local(k) = -ABS( v_comp ) * (                       &
     499              swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                       &
    489500                                     10.0 * ( sk(k,j,i)   - sk(k,j-1,i) )     &
    490501                                   -  5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
     
    495506        ENDIF
    496507       
    497         IF ( i == nxl  .AND.  .NOT. degraded_l )  THEN
     508        IF ( i == i_omp  .AND.  .NOT. degraded_l )  THEN
    498509       
    499510           DO  k = nzb_s_inner(j,i)+1, nzt
    500511              u_comp                 = u(k,j,i) - u_gtrans
    501               swap_flux_x_local(k,j) = u_comp * (                             &
     512              swap_flux_x_local(k,j,tn) = u_comp * (                             &
    502513                                       37.0 * ( sk(k,j,i)   + sk(k,j,i-1) )   &
    503514                                     -  8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
    504515                                     +        ( sk(k,j,i+2) + sk(k,j,i-3) )   &
    505516                                                ) * adv_sca_5
    506               swap_diss_x_local(k,j) = -ABS( u_comp ) * (                     &
     517              swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                     &
    507518                                       10.0 * ( sk(k,j,i)   - sk(k,j,i-1) )   &
    508519                                     -  5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
     
    518529         
    519530          tend(k,j,i) = tend(k,j,i) - (                                      &
    520                           ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
    521                             swap_diss_x_local(k,j) ) * ddx                   &
    522                         + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k)   - &
    523                             swap_diss_y_local(k)   ) * ddy                   &
     531                          ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
     532                            swap_diss_x_local(k,j,tn) ) * ddx                   &
     533                        + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn)   - &
     534                            swap_diss_y_local(k,tn)   ) * ddy                   &
    524535                                      )
    525536
    526           swap_flux_y_local(k)   = flux_n(k)
    527           swap_diss_y_local(k)   = diss_n(k)
    528           swap_flux_x_local(k,j) = flux_r(k)
    529           swap_diss_x_local(k,j) = diss_r(k)
     537          swap_flux_y_local(k,tn)   = flux_n(k)
     538          swap_diss_y_local(k,tn)   = diss_n(k)
     539          swap_flux_x_local(k,j,tn) = flux_r(k)
     540          swap_diss_x_local(k,j,tn) = diss_r(k)
    530541         
    531542       ENDDO
     
    661672! Advection of u-component - Call for grid point i,j
    662673!------------------------------------------------------------------------------!
    663     SUBROUTINE advec_u_ws_ij( i, j )
     674    SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
    664675
    665676       USE arrays_3d
     
    672683       IMPLICIT NONE
    673684
    674        INTEGER ::  i, j, k
     685       INTEGER ::  i, i_omp, j, k, tn
    675686       LOGICAL ::  degraded_l, degraded_s
    676687       REAL    ::  gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp
    677        REAL, DIMENSION(nzb:nzt+1)  ::  flux_t, diss_t, flux_r, diss_r,        &
    678                                        flux_n, diss_n, u_comp
     688       REAL, DIMENSION(nzb:nzt+1) ::  flux_t, diss_t, flux_r, diss_r,        &
     689                                      flux_n, diss_n, u_comp
    679690                                       
    680691       degraded_l = .FALSE.
     
    696707                            -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
    697708                  diss_r(k) = - ABS( u_comp(k) - gu ) * (                     &
    698                               3.0 * ( u(k,j,i+1) - u(k,j,i)   )               & 
     709                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
    699710                            -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
    700711              ENDDO
     
    717728                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    718729                  diss_n(k) = - ABS( v_comp ) * (                             &
    719                               3.0 * ( u(k,j+1,i) - u(k,j,i)     )             & 
     730                              3.0 * ( u(k,j+1,i) - u(k,j,i)     )             &
    720731                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
    721732               ENDDO
     
    735746!--               Compute leftside fluxes for the left boundary of PE domain
    736747                  u_comp(k)     = u(k,j,i) + u(k,j,i-1) - gu
    737                   flux_l_u(k,j) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
    738                   diss_l_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
     748                  flux_l_u(k,j,tn) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
     749                  diss_l_u(k,j,tn) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
    739750                                            u(k,j,i-1), u(k,j,i-1), u_comp(k),&
    740751                                            0.25, ddx )
     
    745756                            -       ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
    746757                  diss_r(k) = - ABS( u_comp(k) -gu ) * (                      &
    747                               3.0 * ( u(k,j,i+1) - u(k,j,i)   )               & 
     758                              3.0 * ( u(k,j,i+1) - u(k,j,i)   )               &
    748759                            -       ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
    749760               ENDDO
     
    757768                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    758769                  diss_n(k) = - ABS( v_comp ) * (                             &
    759                               3.0 * ( u(k,j+1,i) - u(k,j,i)   )               & 
     770                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
    760771                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
    761772               ENDDO
     
    764775               DO  k = nzb_u_inner(j,i)+1, nzt
    765776!               
    766 !--              Compute southside fluxes for the south boundary of PE domain           
     777!--              Compute southside fluxes for the south boundary of PE domain
    767778                  v_comp      = v(k,j,i) + v(k,j,i-1) - gv
    768                   flux_s_u(k) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25
    769                   diss_s_u(k) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i),   &
     779                  flux_s_u(k,tn) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25
     780                  diss_s_u(k,tn) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i),   &
    770781                                          u(k,j-1,i), u(k,j-1,i), v_comp,     &
    771782                                          0.25, ddy )
     
    776787                            -       ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
    777788                  diss_n(k) = - ABS( v_comp ) * (                             &
    778                               3.0 * ( u(k,j+1,i) - u(k,j,i)   )               & 
     789                              3.0 * ( u(k,j+1,i) - u(k,j,i)   )               &
    779790                            -       ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
    780791               ENDDO 
     
    819830       ELSE
    820831!       
    821 !--       Compute the fifth order fluxes for the interior of PE domain.                 
     832!--       Compute the fifth order fluxes for the interior of PE domain.
    822833          DO  k = nzb_u_inner(j,i)+1, nzt
    823834             u_comp(k) = u(k,j,i+1) + u(k,j,i)
     
    844855       ENDIF
    845856!       
    846 !--    Compute left- and southside fluxes for the respective boundary of PE     
     857!--    Compute left- and southside fluxes for the respective boundary of PE
    847858       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
    848859       
    849860          DO  k = nzb_u_inner(j,i)+1, nzt
    850861             v_comp      = v(k,j,i) + v(k,j,i-1) - gv
    851              flux_s_u(k) = v_comp * (                                         &
     862             flux_s_u(k,tn) = v_comp * (                                         &
    852863                           37.0 * ( u(k,j,i) + u(k,j-1,i)   )                 &
    853864                         -  8.0 * ( u(k,j+1,i) + u(k,j-2,i) )                 &
    854865                         +        ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
    855              diss_s_u(k) = - ABS(v_comp) * (                                  &
     866             diss_s_u(k,tn) = - ABS(v_comp) * (                                  &
    856867                           10.0 * ( u(k,j,i) - u(k,j-1,i)   )                 &
    857868                         -  5.0 * ( u(k,j+1,i) - u(k,j-2,i) )                 &
     
    861872       ENDIF
    862873       
    863        IF ( i == nxlu .AND. ( .NOT. degraded_l ) )  THEN
     874       IF ( i == i_omp .AND. ( .NOT. degraded_l ) )  THEN
    864875       
    865876          DO  k = nzb_u_inner(j,i)+1, nzt
    866877             u_comp_l      = u(k,j,i)+u(k,j,i-1)-gu
    867              flux_l_u(k,j) = u_comp_l * (                                     &
     878             flux_l_u(k,j,tn) = u_comp_l * (                                     &
    868879                             37.0 * ( u(k,j,i) + u(k,j,i-1)   )               &
    869880                           -  8.0 * ( u(k,j,i+1) + u(k,j,i-2) )               &
    870881                           +        ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
    871              diss_l_u(k,j) = - ABS(u_comp_l) * (                              &
     882             diss_l_u(k,j,tn) = - ABS(u_comp_l) * (                              &
    872883                             10.0 * ( u(k,j,i) - u(k,j,i-1)   )               &
    873884                           -  5.0 * ( u(k,j,i+1) - u(k,j,i-2) )               &
     
    881892          tend(k,j,i) = tend(k,j,i) - (                                       &
    882893                            ( flux_r(k) + diss_r(k)                           &
    883                           -   flux_l_u(k,j) - diss_l_u(k,j) ) * ddx         &
     894                          -   flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx           &
    884895                          + ( flux_n(k) + diss_n(k)                           &
    885                           -   flux_s_u(k) - diss_s_u(k)     ) * ddy  )
    886 
    887            flux_l_u(k,j) = flux_r(k)
    888            diss_l_u(k,j) = diss_r(k)
    889            flux_s_u(k)   = flux_n(k)
    890            diss_s_u(k)   = diss_n(k)
     896                          -   flux_s_u(k,tn) - diss_s_u(k,tn)     ) * ddy  )
     897
     898           flux_l_u(k,j,tn) = flux_r(k)
     899           diss_l_u(k,j,tn) = diss_r(k)
     900           flux_s_u(k,tn)   = flux_n(k)
     901           diss_s_u(k,tn)   = diss_n(k)
    891902!
    892903!--        Statistical Evaluation of u'u'. The factor has to be applied for
     
    897908             / ( u_comp(k) - gu + 1.0E-20      )                              &
    898909             +   diss_r(k) *                                                  &
    899                  ABS( u_comp(k) - 2.0 * hom(k,1,1,:) )                        & 
     910                 ABS( u_comp(k) - 2.0 * hom(k,1,1,:) )                        &
    900911             / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) )                          &
    901912             *   weight_substep(intermediate_timestep_count) * rmask(j,i,:)
     
    9981009
    9991010
    1000 !------------------------------------------------------------------------------!
     1011!-----------------------------------------------------------------------------!
    10011012! Advection of v-component - Call for grid point i,j
    1002 !------------------------------------------------------------------------------!
    1003    SUBROUTINE advec_v_ws_ij( i, j )
     1013!-----------------------------------------------------------------------------!
     1014   SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
    10041015
    10051016       USE arrays_3d
     
    10121023       IMPLICIT NONE
    10131024
    1014        INTEGER  ::  i, j, k
     1025       INTEGER  ::  i, i_omp, j, k, tn
    10151026       LOGICAL  ::  degraded_l, degraded_s
    10161027       REAL     ::  gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp
     
    10361047                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    10371048                   diss_r(k) = - ABS( u_comp ) * (                            &
    1038                                3.0 * ( v(k,j,i+1) - v(k,j,i)   )              & 
     1049                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
    10391050                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
    10401051                ENDDO
     
    10561067                             -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
    10571068                   diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
    1058                                3.0 * ( v(k,j+1,i) - v(k,j,i)   )              & 
     1069                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
    10591070                             -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
    10601071                ENDDO
     
    10771088                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    10781089                   diss_r(k) = - ABS( u_comp ) * (                            &
    1079                                3.0 * ( w(k,j,i+1) - w(k,j,i)   )              & 
     1090                               3.0 * ( w(k,j,i+1) - w(k,j,i)   )              &
    10801091                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
    10811092                ENDDO
     
    10841095                DO  k = nzb_v_inner(j,i)+1, nzt
    10851096                   u_comp        = u(k,j-1,i) + u(k,j,i) - gu
    1086                    flux_l_v(k,j) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
    1087                    diss_l_v(k,j) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),&
     1097                   flux_l_v(k,j,tn) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
     1098                   diss_l_v(k,j,tn) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),&
    10881099                                             v(k,j,i-1), v(k,j,i-1), u_comp,  &
    10891100                                             0.25, ddx )
     
    10941105                             -       ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
    10951106                   diss_r(k) = - ABS( u_comp ) * (                            &
    1096                                3.0 * ( v(k,j,i+1) - v(k,j,i)   )              & 
     1107                               3.0 * ( v(k,j,i+1) - v(k,j,i)   )              &
    10971108                             -       ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
    10981109                ENDDO
     
    11021113                DO  k = nzb_v_inner(j,i)+1, nzt
    11031114                   v_comp(k)   = v(k,j,i) + v(k,j-1,i) - gv
    1104                    flux_s_v(k) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
    1105                    diss_s_v(k) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i),  &
     1115                   flux_s_v(k,tn) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
     1116                   diss_s_v(k,tn) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i),  &
    11061117                                           v(k,j-1,i), v(k,j-1,i), v_comp(k), &
    11071118                                           0.25, ddy )
     
    11121123                             -       ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
    11131124                   diss_n(k) = - ABS( v_comp(k) - gv ) * (                    &
    1114                                3.0 * ( v(k,j+1,i) - v(k,j,i)   )              & 
     1125                               3.0 * ( v(k,j+1,i) - v(k,j,i)   )              &
    11151126                             -       ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
    11161127                ENDDO
     
    11811192!       
    11821193!--    Compute left- and southside fluxes for the respective boundary       
    1183        IF ( i == nxl .AND. ( .NOT. degraded_l ) )  THEN
     1194       IF ( i == i_omp  .AND.  .NOT. degraded_l )  THEN
     1195
    11841196          DO  k = nzb_v_inner(j,i)+1, nzt
    11851197             u_comp        = u(k,j-1,i) + u(k,j,i) - gu
    1186              flux_l_v(k,j) = u_comp * (                                       &
     1198             flux_l_v(k,j,tn) = u_comp * (                                       &
    11871199                             37.0 * ( v(k,j,i) + v(k,j,i-1)   )               &
    11881200                           -  8.0 * ( v(k,j,i+1) + v(k,j,i-2) )               &
    11891201                           +        ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
    1190              diss_l_v(k,j) = - ABS( u_comp ) * (                              &
     1202             diss_l_v(k,j,tn) = - ABS( u_comp ) * (                              &
    11911203                             10.0 * ( v(k,j,i) - v(k,j,i-1)   )               &
    11921204                           -  5.0 * ( v(k,j,i+1) - v(k,j,i-2) )               &
     
    11961208       ENDIF
    11971209       
    1198        IF ( j == nysv .AND. ( .NOT. degraded_s ) )  THEN
     1210       IF ( j == nysv  .AND.  .NOT. degraded_s )  THEN
    11991211       
    12001212          DO  k = nzb_v_inner(j,i)+1, nzt
    12011213             v_comp_l    = v(k,j,i) + v(k,j-1,i) - gv
    1202              flux_s_v(k) = v_comp_l * (                                       &
     1214             flux_s_v(k,tn) = v_comp_l * (                                       &
    12031215                           37.0 * ( v(k,j,i) + v(k,j-1,i)   )                 &
    12041216                         -  8.0 * ( v(k,j+1,i) + v(k,j-2,i) )                 &
    12051217                         +        ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
    1206              diss_s_v(k) = - ABS( v_comp_l ) * (                              &
     1218             diss_s_v(k,tn) = - ABS( v_comp_l ) * (                              &
    12071219                           10.0 * ( v(k,j,i) - v(k,j-1,i)   )                 &
    12081220                         -  5.0 * ( v(k,j+1,i) - v(k,j-2,i) )                 &
     
    12161228          tend(k,j,i) = tend(k,j,i) - (                                       &
    12171229                         ( flux_r(k) + diss_r(k)                              &
    1218                        -   flux_l_v(k,j) - diss_l_v(k,j)   ) * ddx            &
     1230                       -   flux_l_v(k,j,tn) - diss_l_v(k,j,tn)   ) * ddx            &
    12191231                       + ( flux_n(k) + diss_n(k)                              &
    1220                        -   flux_s_v(k) - diss_s_v(k)       ) * ddy     )
     1232                       -   flux_s_v(k,tn) - diss_s_v(k,tn)       ) * ddy     )
    12211233       
    1222            flux_l_v(k,j) = flux_r(k)
    1223            diss_l_v(k,j) = diss_r(k)
    1224            flux_s_v(k)   = flux_n(k)
    1225            diss_s_v(k)   = diss_n(k)
     1234           flux_l_v(k,j,tn) = flux_r(k)
     1235           diss_l_v(k,j,tn) = diss_r(k)
     1236           flux_s_v(k,tn)   = flux_n(k)
     1237           diss_s_v(k,tn)   = diss_n(k)
    12261238
    12271239!
     
    13331345! Advection of w-component - Call for grid point i,j
    13341346!------------------------------------------------------------------------------!
    1335     SUBROUTINE advec_w_ws_ij( i, j )
     1347    SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
    13361348
    13371349       USE arrays_3d
     
    13441356       IMPLICIT NONE
    13451357
    1346        INTEGER ::  i, j, k
     1358       INTEGER ::  i, i_omp, j, k, tn
    13471359       LOGICAL ::  degraded_l, degraded_s
    13481360       REAL    ::  gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp
     
    13551367       gu = 2.0 * u_gtrans
    13561368       gv = 2.0 * v_gtrans
    1357        
     1369
     1370             
    13581371       IF ( boundary_flags(j,i) /= 0 )  THEN
    13591372!       
     
    14251438                 
    14261439                  u_comp        = u(k+1,j,i) + u(k,j,i) - gu
    1427                   flux_l_w(k,j) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25
    1428                   diss_l_w(k,j) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
     1440                  flux_l_w(k,j,tn) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25
     1441                  diss_l_w(k,j,tn) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
    14291442                                            w(k,j,i-1), w(k,j,i-1), u_comp,   &
    14301443                                            0.25, ddx )
     
    14551468!--              Compute southside fluxes for the south boundary of PE domain           
    14561469                  v_comp      = v(k+1,j,i) + v(k,j,i) - gv
    1457                   flux_s_w(k) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25
    1458                   diss_s_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i),   &
     1470                  flux_s_w(k,tn) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25
     1471                  diss_s_w(k,tn) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i),   &
    14591472                                          w(k,j-1,i), w(k,j-1,i), v_comp,     &
    14601473                                          0.25, ddy )                 
     
    15261539!       
    15271540!--    Compute left- and southside fluxes for the respective boundary     
    1528        IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
     1541       IF ( j == nys  .AND.  .NOT. degraded_s )  THEN
    15291542       
    15301543          DO  k = nzb_w_inner(j,i)+1, nzt
    15311544             v_comp      = v(k+1,j,i) + v(k,j,i) - gv
    1532              flux_s_w(k) = v_comp * (                                         &
     1545             flux_s_w(k,tn) = v_comp * (                                         &
    15331546                           37.0 * ( w(k,j,i) + w(k,j-1,i)   )                 &
    15341547                         -  8.0 * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
    15351548                         +        ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
    1536              diss_s_w(k) = - ABS( v_comp ) * (                                &
     1549             diss_s_w(k,tn) = - ABS( v_comp ) * (                                &
    15371550                           10.0 * ( w(k,j,i) - w(k,j-1,i)   )                 &
    15381551                         -  5.0 * ( w(k,j+1,i) - w(k,j-2,i) )                 &
     
    15421555       ENDIF
    15431556       
    1544        IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN
    1545        
     1557       IF ( i == i_omp  .AND.   .NOT. degraded_l ) THEN
     1558          
    15461559          DO  k = nzb_w_inner(j,i)+1, nzt
    15471560            u_comp        = u(k+1,j,i) + u(k,j,i) - gu
    1548             flux_l_w(k,j) = u_comp * (                                        &
     1561            flux_l_w(k,j,tn) = u_comp * (                                        &
    15491562                            37.0 * ( w(k,j,i) + w(k,j,i-1)   )                &
    15501563                          -  8.0 * ( w(k,j,i+1) + w(k,j,i-2) )                &
    15511564                          +        ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
    1552             diss_l_w(k,j) = - ABS( u_comp ) * (                               &
     1565            diss_l_w(k,j,tn) = - ABS( u_comp ) * (                               &
    15531566                            10.0 * ( w(k,j,i) - w(k,j,i-1)   )                &
    15541567                          -  5.0 * ( w(k,j,i+1) - w(k,j,i-2) )                &
     
    15571570         
    15581571       ENDIF
     1572
    15591573!       
    15601574!--    Now compute the tendency terms for the horizontal parts.         
     
    15621576          tend(k,j,i) = tend(k,j,i) - (                                       &
    15631577                      ( flux_r(k) + diss_r(k)                                 &
    1564                     -   flux_l_w(k,j) - diss_l_w(k,j)   ) * ddx               &
     1578                    -   flux_l_w(k,j,tn) - diss_l_w(k,j,tn)   ) * ddx               &
    15651579                    + ( flux_n(k) + diss_n(k)                                 &
    1566                     -   flux_s_w(k) - diss_s_w(k)       ) * ddy  )
    1567 
    1568           flux_l_w(k,j) = flux_r(k)
    1569           diss_l_w(k,j) = diss_r(k)
    1570           flux_s_w(k) = flux_n(k)
    1571           diss_s_w(k) = diss_n(k)
     1580                    -   flux_s_w(k,tn) - diss_s_w(k,tn)       ) * ddy  )
     1581
     1582          flux_l_w(k,j,tn) = flux_r(k)
     1583          diss_l_w(k,j,tn) = diss_r(k)
     1584          flux_s_w(k,tn) = flux_n(k)
     1585          diss_s_w(k,tn) = diss_n(k)
    15721586       ENDDO
    15731587
Note: See TracChangeset for help on using the changeset viewer.