Ignore:
Timestamp:
Jul 9, 2019 6:04:41 PM (2 years ago)
Author:
suehring
Message:

Implementation of a monotonic flux limiter for vertical advection term in Wicker-Skamarock scheme. The flux limiter is currently only applied for passive scalars (passive scalar, chemical species, aerosols) within the region up to the highest topography, in order to avoid the built-up of large concentrations within poorly resolved cavities in urban environments. To enable the limiter monotonic_limiter_z = .T. must be set. Note, the limiter is currently only implemented for the cache-optimized version of advec_ws. Further changes in offline nesting: Set boundary condition for w at nzt+1 at all lateral boundaries (even though these won't enter the numerical solution), in order to avoid high vertical velocities in the run-control file which might built-up due to the mass-conservation; bugfix in offline nesting for chemical species

File:
1 edited

Legend:

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

    r3873 r4079  
    2525! -----------------
    2626! $Id$
     27! Implementation of a flux limiter according to Skamarock (2006) for the
     28! vertical scalar advection. Please note, this is only implemented for the
     29! cache-optimized version at the moment. Implementation for the vector-
     30! optimized version will follow after critical issues concerning
     31! vectorization are fixed.
     32!
     33! 3873 2019-04-08 15:44:30Z knoop
    2734! Moved ocean_mode specific code to ocean_mod
    2835!
     
    273280!> the divergence is not sufficiently reduced, resulting in erroneous fluxes
    274281!> and could lead to numerical instabilities.
     282!>
     283!> @todo Implement monotonic flux limiter also for vector version.
    275284!------------------------------------------------------------------------------!
    276285 MODULE advec_ws
     
    296305               bc_dirichlet_s, bc_radiation_l, bc_radiation_n,                 &
    297306               bc_radiation_r, bc_radiation_s, intermediate_timestep_count,    &
    298                u_gtrans, v_gtrans
     307               u_gtrans, v_gtrans, dt_3d
    299308
    300309    USE indices,                                                               &
     
    11161125    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local,            &
    11171126                              swap_diss_y_local, swap_flux_x_local,            &
    1118                               swap_diss_x_local, i_omp, tn )
     1127                              swap_diss_x_local, i_omp, tn, flux_limitation )
    11191128
    11201129
     
    11261135       INTEGER(iwp) ::  k         !< grid index along z-direction
    11271136       INTEGER(iwp) ::  k_mm      !< k-2 index in disretization, can be modified to avoid segmentation faults
     1137       INTEGER(iwp) ::  k_mmm     !< k-3 index in disretization, can be modified to avoid segmentation faults
    11281138       INTEGER(iwp) ::  k_pp      !< k+2 index in disretization, can be modified to avoid segmentation faults
    11291139       INTEGER(iwp) ::  k_ppp     !< k+3 index in disretization, can be modified to avoid segmentation faults
    11301140       INTEGER(iwp) ::  nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 
    11311141       INTEGER(iwp) ::  tn        !< number of OpenMP thread
    1132        
    1133        REAL(wp)     ::  ibit0  !< flag indicating 1st-order scheme along x-direction
    1134        REAL(wp)     ::  ibit1  !< flag indicating 3rd-order scheme along x-direction
    1135        REAL(wp)     ::  ibit2  !< flag indicating 5th-order scheme along x-direction
    1136        REAL(wp)     ::  ibit3  !< flag indicating 1st-order scheme along y-direction
    1137        REAL(wp)     ::  ibit4  !< flag indicating 3rd-order scheme along y-direction
    1138        REAL(wp)     ::  ibit5  !< flag indicating 5th-order scheme along y-direction
    1139        REAL(wp)     ::  ibit6  !< flag indicating 1st-order scheme along z-direction
    1140        REAL(wp)     ::  ibit7  !< flag indicating 3rd-order scheme along z-direction
    1141        REAL(wp)     ::  ibit8  !< flag indicating 5th-order scheme along z-direction
    1142        REAL(wp)     ::  diss_d !< artificial dissipation term at grid box bottom
    1143        REAL(wp)     ::  div    !< diverence on scalar grid
    1144        REAL(wp)     ::  flux_d !< 6th-order flux at grid box bottom
    1145        REAL(wp)     ::  u_comp !< advection velocity along x-direction
    1146        REAL(wp)     ::  v_comp !< advection velocity along y-direction
     1142
     1143       LOGICAL, OPTIONAL ::  flux_limitation !< flag indicating flux limitation of the vertical advection
     1144       LOGICAL           ::  limiter         !< control flag indicating the application of flux limitation
     1145
     1146       REAL(wp) ::  diss_d        !< artificial dissipation term at grid box bottom
     1147       REAL(wp) ::  div           !< velocity diverence on scalar grid
     1148       REAL(wp) ::  div_in        !< vertical flux divergence of ingoing fluxes
     1149       REAL(wp) ::  div_out       !< vertical flux divergence of outgoing fluxes     
     1150       REAL(wp) ::  f_corr_t      !< correction flux at grid-cell top, i.e. the difference between high and low-order flux
     1151       REAL(wp) ::  f_corr_d      !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux
     1152       REAL(wp) ::  f_corr_t_in   !< correction flux of ingoing flux part at grid-cell top
     1153       REAL(wp) ::  f_corr_d_in   !< correction flux of ingoing flux part at grid-cell bottom
     1154       REAL(wp) ::  f_corr_t_out  !< correction flux of outgoing flux part at grid-cell top
     1155       REAL(wp) ::  f_corr_d_out  !< correction flux of outgoing flux part at grid-cell bottom
     1156       REAL(wp) ::  fac_correction!< factor to limit the in- and outgoing fluxes
     1157       REAL(wp) ::  flux_d        !< 6th-order flux at grid box bottom
     1158       REAL(wp) ::  ibit0         !< flag indicating 1st-order scheme along x-direction
     1159       REAL(wp) ::  ibit1         !< flag indicating 3rd-order scheme along x-direction
     1160       REAL(wp) ::  ibit2         !< flag indicating 5th-order scheme along x-direction
     1161       REAL(wp) ::  ibit3         !< flag indicating 1st-order scheme along y-direction
     1162       REAL(wp) ::  ibit4         !< flag indicating 3rd-order scheme along y-direction
     1163       REAL(wp) ::  ibit5         !< flag indicating 5th-order scheme along y-direction
     1164       REAL(wp) ::  ibit6         !< flag indicating 1st-order scheme along z-direction
     1165       REAL(wp) ::  ibit7         !< flag indicating 3rd-order scheme along z-direction
     1166       REAL(wp) ::  ibit8         !< flag indicating 5th-order scheme along z-direction
     1167       REAL(wp) ::  max_val       !< maximum value of the quanitity along the numerical stencil (in vertical direction)
     1168       REAL(wp) ::  min_val       !< maximum value of the quanitity along the numerical stencil (in vertical direction)
     1169       REAL(wp) ::  mon           !< monotone solution of the advection equation using 1st-order fluxes
     1170       REAL(wp) ::  u_comp        !< advection velocity along x-direction
     1171       REAL(wp) ::  v_comp        !< advection velocity along y-direction
    11471172!       
    11481173!--    sk is an array from parameter list. It should not be a pointer, because
     
    11531178       REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<  advected scalar
    11541179
    1155        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n !< discretized artificial dissipation at northward-side of the grid box
    1156        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r !< discretized artificial dissipation at rightward-side of the grid box
    1157        REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_t !< discretized artificial dissipation at top of the grid box
    1158        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_n !< discretized 6th-order flux at northward-side of the grid box
    1159        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r !< discretized 6th-order flux at rightward-side of the grid box
    1160        REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t !< discretized 6th-order flux at top of the grid box
     1180       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_n     !< discretized artificial dissipation at northward-side of the grid box
     1181       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_r     !< discretized artificial dissipation at rightward-side of the grid box
     1182       REAL(wp), DIMENSION(nzb:nzt+1)         ::  diss_t     !< discretized artificial dissipation at top of the grid box
     1183       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_n     !< discretized 6th-order flux at northward-side of the grid box
     1184       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_r     !< discretized 6th-order flux at rightward-side of the grid box
     1185       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t     !< discretized 6th-order flux at top of the grid box
     1186       REAL(wp), DIMENSION(nzb:nzt+1)         ::  flux_t_1st !< discretized 1st-order flux at top of the grid box
    11611187       
    11621188       REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) ::  swap_diss_y_local !< discretized artificial dissipation at southward-side of the grid box
     
    11781204          nzb_max_l = nzb_max
    11791205       END IF
     1206!
     1207!--    Set control flag for flux limiter
     1208       limiter = .FALSE.
     1209       IF ( PRESENT( flux_limitation) )  limiter = flux_limitation
    11801210!
    11811211!--    Compute southside fluxes of the respective PE bounds.
     
    12311261                                  +           ( sk(k,j+2,i) + sk(k,j-3,i) )   &
    12321262                                                ) * adv_sca_5
    1233               swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                    &
     1263             swap_diss_y_local(k,tn) = -ABS( v_comp ) * (                     &
    12341264                                    10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )   &
    12351265                                  -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )   &
    12361266                                  +             sk(k,j+2,i) - sk(k,j-3,i)     &
    1237                                                          ) * adv_sca_5
     1267                                                        ) * adv_sca_5
    12381268
    12391269          ENDDO
     
    12661296                                                  )
    12671297
    1268               swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                  &
     1298             swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * (                   &
    12691299                                               ( 10.0_wp * ibit2 * adv_sca_5  &
    12701300                                            +     3.0_wp * ibit1 * adv_sca_3  &
     
    12791309                                               ) *                            &
    12801310                                            ( sk(k,j,i+2) - sk(k,j,i-3)    )  &
    1281                                                            )
     1311                                                          )
    12821312
    12831313          ENDDO
     
    15351565                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    15361566                                                         )
    1537        ENDDO
     1567       ENDDO
     1568
     1569       IF ( limiter )  THEN
     1570!
     1571!--       Compute monotone first-order fluxes which are required for mononte
     1572!--       flux limitation.
     1573          flux_t_1st(nzb) = 0.0_wp
     1574          DO  k = nzb+1, nzb_max_l
     1575             flux_t_1st(k) = ( w(k,j,i)   * ( sk(k+1,j,i)  + sk(k,j,i) )       &
     1576                       -  ABS( w(k,j,i) ) * ( sk(k+1,j,i)  - sk(k,j,i) ) )     &
     1577                           * rho_air_zw(k) * adv_sca_1
     1578!
     1579!--          In flux limitation the total flux will be corrected. For the sake
     1580!--          of cleariness the higher-order advective and disspative fluxes
     1581!--          will be merged onto flux_t.
     1582             flux_t(k) = flux_t(k) + diss_t(k)
     1583             diss_t(k) = 0.0_wp
     1584          ENDDO
     1585!
     1586!--       Flux limitation of vertical fluxes according to Skamarock (2006).
     1587!--       Please note, as flux limitation implies linear dependencies of fluxes,
     1588!--       flux limitation is only made for the vertical advection term. Limitation
     1589!--       of the horizontal terms cannot be parallelized.
     1590!--       Due to the linear dependency, the following loop will not be vectorized.
     1591!--       Further, note that the flux limiter is only applied within the urban
     1592!--       layer, i.e up to the topography top. 
     1593          DO  k = nzb+1, nzb_max_l
     1594!
     1595!--          Compute one-dimensional divergence along the vertical direction,
     1596!--          which is used to correct the advection discretization. This is
     1597!--          necessary as in one-dimensional space the advection velocity
     1598!--          should actually be constant.
     1599             div = ( w(k,j,i)   * rho_air_zw(k)                                &
     1600                   - w(k-1,j,i) * rho_air_zw(k-1)                              &     
     1601                   ) * drho_air(k) * ddzw(k)
     1602!
     1603!--          Compute monotone solution of the advection equation from
     1604!--          1st-order fluxes. Please note, the advection equation is corrected
     1605!--          by the divergence term (in 1D the advective flow should be divergence
     1606!--          free). Moreover, please note, as time-increment the full timestep
     1607!--          is used, even though a Runge-Kutta scheme will be used. However,
     1608!--          the length of the actual time increment is not important at all
     1609!--          since it cancels out later when the fluxes are limited.   
     1610             mon = sk(k,j,i) + ( - ( flux_t_1st(k) - flux_t_1st(k-1) )         &
     1611                             * drho_air(k) * ddzw(k)                           &
     1612                             + div * sk(k,j,i)                                 &
     1613                               ) * dt_3d
     1614!
     1615!--          Determine minimum and maximum values along the numerical stencil.
     1616             k_mmm = MAX( k - 3, nzb + 1 )
     1617             k_ppp = MIN( k + 3, nzt + 1 )
     1618
     1619             min_val = MINVAL( sk(k_mmm:k_ppp,j,i) )
     1620             max_val = MAXVAL( sk(k_mmm:k_ppp,j,i) )
     1621!
     1622!--          Compute difference between high- and low-order fluxes, which may
     1623!--          act as correction fluxes
     1624             f_corr_t = flux_t(k)   - flux_t_1st(k)
     1625             f_corr_d = flux_t(k-1) - flux_t_1st(k-1)
     1626!
     1627!--          Determine outgoing fluxes, i.e. the part of the fluxes which can
     1628!--          decrease the value within the grid box
     1629             f_corr_t_out = MAX( 0.0_wp, f_corr_t )
     1630             f_corr_d_out = MIN( 0.0_wp, f_corr_d )
     1631!
     1632!--          Determine ingoing fluxes, i.e. the part of the fluxes which can
     1633!--          increase the value within the grid box
     1634             f_corr_t_in = MIN( 0.0_wp, f_corr_t)
     1635             f_corr_d_in = MAX( 0.0_wp, f_corr_d)
     1636!
     1637!--          Compute divergence of outgoing correction fluxes
     1638             div_out = - ( f_corr_t_out - f_corr_d_out ) * drho_air(k)         &
     1639                                                         * ddzw(k) * dt_3d
     1640!
     1641!--          Compute divergence of ingoing correction fluxes
     1642             div_in = - ( f_corr_t_in - f_corr_d_in )    * drho_air(k)         &
     1643                                                         * ddzw(k) * dt_3d
     1644!
     1645!--          Check if outgoing fluxes can lead to undershoots, i.e. values smaller
     1646!--          than the minimum value within the numerical stencil. If so, limit
     1647!--          them.
     1648             IF ( mon - min_val < - div_out  .AND.  ABS( div_out ) > 0.0_wp )  &
     1649             THEN
     1650                fac_correction = ( mon - min_val ) / ( - div_out )
     1651                f_corr_t_out = f_corr_t_out * fac_correction
     1652                f_corr_d_out = f_corr_d_out * fac_correction
     1653             ENDIF
     1654!
     1655!--          Check if ingoing fluxes can lead to overshoots, i.e. values larger
     1656!--          than the maximum value within the numerical stencil. If so, limit
     1657!--          them.
     1658             IF ( mon - max_val > - div_in  .AND.  ABS( div_in ) > 0.0_wp )    &
     1659             THEN
     1660                fac_correction = ( mon - max_val ) / ( - div_in )
     1661                f_corr_t_in = f_corr_t_in * fac_correction
     1662                f_corr_d_in = f_corr_d_in * fac_correction
     1663             ENDIF
     1664!
     1665!--          Finally add the limited fluxes to the original ones. If no
     1666!--          flux limitation was done, the fluxes equal the original ones.
     1667             flux_t(k)   = flux_t_1st(k)   + f_corr_t_out + f_corr_t_in
     1668             flux_t(k-1) = flux_t_1st(k-1) + f_corr_d_out + f_corr_d_in
     1669          ENDDO
     1670       ENDIF
    15381671
    15391672       DO  k = nzb+1, nzb_max_l
Note: See TracChangeset for help on using the changeset viewer.