Changeset 736 for palm/trunk/SOURCE


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

Bugfix in ws-scheme concerning OpenMP paralellization.

Location:
palm/trunk/SOURCE
Files:
3 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
  • palm/trunk/SOURCE/modules.f90

    r723 r736  
    55! Current revisions:
    66! -----------------
    7 !
     7! Dimension of fluxes needed for WS-scheme increased.
    88!
    99! Former revisions:
     
    246246          km_damp_x, km_damp_y, lad, l_grid, pt_init, q_init, rdf, sa_init,    &
    247247          ug, u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs,  &
    248           zu, zw, flux_s_u, flux_s_v, flux_s_w, diss_s_u, diss_s_v, diss_s_w,  &
    249           flux_s_pt, diss_s_pt, flux_s_e, diss_s_e, flux_s_q, diss_s_q,        &
    250           flux_s_sa, diss_s_sa
     248          zu, zw
    251249
    252250    REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
    253251          c_u, c_v, c_w, dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg,                  &
    254252          mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, ts, us, z0,    &
    255           flux_l_u, flux_l_v, flux_l_w, diss_l_u, diss_l_v, diss_l_w,          &
    256           flux_l_pt, diss_l_pt, flux_l_e, diss_l_e, flux_l_q, diss_l_q,        &
    257           flux_l_sa, diss_l_sa, total_2d_o, total_2d_a
     253          flux_s_pt, diss_s_pt, flux_s_e, diss_s_e, flux_s_q, diss_s_q,        &
     254          flux_s_sa, diss_s_sa, flux_s_u, flux_s_v, flux_s_w, diss_s_u,        &
     255          diss_s_v, diss_s_w, total_2d_o, total_2d_a
    258256
    259257    REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
     
    270268          canopy_heat_flux, cdc, d, diss, lad_s, lad_u, lad_v, lad_w, lai,     &
    271269          l_wall, p_loc, sec, sls, tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l,    &
    272           v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
     270          v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s, flux_l_pt,          &
     271          diss_l_pt, flux_l_e, diss_l_e, flux_l_q, diss_l_q, flux_l_sa,        &
     272          diss_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_l_u, diss_l_v, diss_l_w
    273273
    274274    REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
  • palm/trunk/SOURCE/prognostic_equations.f90

    r710 r736  
    44! Current revisions:
    55! -----------------
    6 !
     6! Bugfix: determination of first thread index i for WS-scheme
    77!
    88! Former revisions:
     
    147147
    148148    CHARACTER (LEN=9) ::  time_to_string
    149     INTEGER ::  i, j, k
     149    INTEGER ::  i, i_omp_start, j, k, tn = 0
    150150    REAL    ::  sat, sbt
    151151
     
    172172!
    173173!-- u-tendency terms with no communication
     174    i_omp_start = nxlu
    174175    DO  i = nxlu, nxr
    175176       DO  j = nys, nyn
     
    179180             tend(:,j,i) = 0.0
    180181             IF ( ws_scheme_mom )  THEN
    181                  CALL advec_u_ws( i, j )
     182                 CALL advec_u_ws( i, j, i_omp_start, tn )
    182183             ELSE
    183184                 CALL advec_u_pw( i, j )
     
    257258!
    258259!-- v-tendency terms with no communication
     260    i_omp_start = nxl
    259261    DO  i = nxl, nxr
    260262       DO  j = nysv, nyn
     
    264266             tend(:,j,i) = 0.0
    265267             IF ( ws_scheme_mom )  THEN
    266                  CALL advec_v_ws( i, j )
     268                 CALL advec_v_ws( i, j, i_omp_start, tn )
    267269             ELSE
    268270                 CALL advec_v_pw( i, j )
     
    348350             tend(:,j,i) = 0.0
    349351             IF ( ws_scheme_mom )  THEN
    350                  CALL advec_w_ws( i, j )
     352                 CALL advec_w_ws( i, j, i_omp_start, tn )
    351353             ELSE 
    352354                 CALL advec_w_pw( i, j )
     
    453455                IF ( ws_scheme_sca )  THEN
    454456                   CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
    455                                  diss_s_pt, flux_l_pt, diss_l_pt )
     457                                 diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
    456458                ELSE
    457459                    CALL advec_s_pw( i, j, pt )
     
    574576                   IF ( ws_scheme_sca )  THEN
    575577                       CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
    576                                     diss_s_sa, flux_l_sa, diss_l_sa )
     578                                    diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn )
    577579                   ELSE
    578580                       CALL advec_s_pw( i, j, sa )
     
    673675                   IF ( ws_scheme_sca )  THEN
    674676                       CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
    675                                    diss_s_q, flux_l_q, diss_l_q )
     677                                   diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
    676678                   ELSE
    677679                       CALL advec_s_pw( i, j, q )
     
    813815                      IF ( ws_scheme_sca )  THEN
    814816                          CALL advec_s_ws( i, j, e, 'e', flux_s_e, &
    815                                       diss_s_e, flux_l_e, diss_l_e )
     817                                      diss_s_e, flux_l_e, diss_l_e, i_omp_start, tn )
    816818                      ELSE
    817819                          CALL advec_s_pw( i, j, e )
     
    916918
    917919    CHARACTER (LEN=9) ::  time_to_string
    918     INTEGER ::  i, j, k
     920    INTEGER ::  i, i_omp_start, j, k, omp_get_thread_num, tn = 0
     921    LOGICAL ::  loop_start
    919922
    920923
     
    934937         intermediate_timestep_count == 1 )  CALL ws_statistics
    935938
    936 
    937939!
    938940!-- Loop over all prognostic equations
    939 !$OMP PARALLEL private (i,j,k)
     941!$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
     942
     943!$  tn = omp_get_thread_num()
     944    loop_start = .TRUE.
    940945!$OMP DO
    941946    DO  i = nxl, nxr
     947
     948!
     949!--    Store the first loop index. It differs for each thread and is required
     950!--    later in advec_ws
     951       IF ( loop_start )  THEN
     952          loop_start  = .FALSE.
     953          i_omp_start = i
     954       ENDIF
     955 
    942956       DO  j = nys, nyn
    943957!
     
    948962             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
    949963                IF ( ws_scheme_mom )  THEN
    950     !               CALL local_diss( i, j, u)    ! dissipation control
    951                    CALL advec_u_ws( i, j )
     964                   IF ( outflow_l .AND. i_omp_start == nxl )  THEN
     965!                     CALL local_diss( i, j, u)    ! dissipation control
     966                      CALL advec_u_ws( i, j, i_omp_start + 1, tn )
     967                   ELSE
     968                      CALL advec_u_ws( i, j, i_omp_start, tn )
     969                   ENDIF
    952970                ELSE
    953971                   CALL advec_u_pw( i, j )
    954972                ENDIF
    955             ELSE
     973             ELSE
    956974                CALL advec_u_up( i, j )
    957975             ENDIF
     
    10171035                IF ( ws_scheme_mom )  THEN
    10181036                 !   CALL local_diss( i, j, v)
    1019                     CALL advec_v_ws( i, j )
     1037                    CALL advec_v_ws( i, j, i_omp_start, tn )
    10201038                ELSE
    10211039                    CALL advec_v_pw( i, j )
     
    10811099             IF ( ws_scheme_mom )  THEN
    10821100             !   CALL local_diss( i, j, w)
    1083                 CALL advec_w_ws( i, j )
     1101                CALL advec_w_ws( i, j, i_omp_start, tn )
    10841102             ELSE
    10851103                CALL advec_w_pw( i, j )
     
    11451163       !            CALL local_diss( i, j, pt )
    11461164                   CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
    1147                              diss_s_pt, flux_l_pt, diss_l_pt )
     1165                             diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
    11481166                ELSE
    11491167                   CALL advec_s_pw( i, j, pt )
     
    12271245            !        CALL local_diss( i, j, sa )
    12281246                    CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
    1229                                 diss_s_sa, flux_l_sa, diss_l_sa  )
     1247                                diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn  )
    12301248                ELSE
    12311249                    CALL advec_s_pw( i, j, sa )
     
    12851303          !         CALL local_diss( i, j, q )
    12861304                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
    1287                                 diss_s_q, flux_l_q, diss_l_q )
     1305                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
    12881306                ELSE
    12891307                   CALL advec_s_pw( i, j, q )
     
    13611379                 !    CALL local_diss( i, j, e )
    13621380                     CALL advec_s_ws( i, j, e, 'e', flux_s_e, &
    1363                                 diss_s_e, flux_l_e, diss_l_e )
     1381                                diss_s_e, flux_l_e, diss_l_e , i_omp_start, tn )
    13641382                 ELSE
    13651383                     CALL advec_s_pw( i, j, e )
Note: See TracChangeset for help on using the changeset viewer.