Ignore:
Timestamp:
Apr 11, 2017 11:37:51 AM (4 years ago)
Author:
suehring
Message:

Bugfixes in radiation_model; monotonic limiter in 5th order advection scheme removed

File:
1 edited

Legend:

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

    r2119 r2200  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! monotonic_adjustment removed
    2323!
    2424! Former revisions:
     
    245245       USE control_parameters,                                                 &
    246246           ONLY:  cloud_physics, humidity, loop_optimization,                  &
    247                   monotonic_adjustment, passive_scalar, microphysics_seifert,  &
     247                  passive_scalar, microphysics_seifert,                        &
    248248                  ocean, ws_scheme_mom, ws_scheme_sca
    249249
     
    258258           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
    259259                  sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l,            &
    260                   sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsss_ws_l,               &
     260                  sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsss_ws_l,            &
    261261                  sums_wsus_ws_l, sums_wsvs_ws_l 
    262262
     
    418418   
    419419
    420        IF ( scalar_advec == 'ws-scheme' .OR.                                   &
    421             scalar_advec == 'ws-scheme-mono' )  THEN
     420       IF ( scalar_advec == 'ws-scheme' )  THEN
    422421!
    423422!--       Set flags to steer the degradation of the advection scheme in advec_ws
     
    803802!
    804803!--    Exchange 3D integer wall_flags.
    805        IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme'  &
    806        .OR. scalar_advec == 'ws-scheme-mono' )  THEN 
     804       IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme'     &
     805         )  THEN 
    807806!
    808807!--       Exchange ghost points for advection flags
     
    901900
    902901       USE control_parameters,                                                 &
    903            ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans, &
    904                   v_gtrans
     902           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    905903
    906904       USE grid_variables,                                                     &
     
    946944       REAL(wp)     ::  div    !<
    947945       REAL(wp)     ::  flux_d !<
    948        REAL(wp)     ::  fd_1   !<
    949        REAL(wp)     ::  fl_1   !<
    950        REAL(wp)     ::  fn_1   !<
    951        REAL(wp)     ::  fr_1   !<
    952        REAL(wp)     ::  fs_1   !<
    953        REAL(wp)     ::  ft_1   !<
    954        REAL(wp)     ::  phi_d  !<
    955        REAL(wp)     ::  phi_l  !<
    956        REAL(wp)     ::  phi_n  !<
    957        REAL(wp)     ::  phi_r  !<
    958        REAL(wp)     ::  phi_s  !<
    959        REAL(wp)     ::  phi_t  !<
    960        REAL(wp)     ::  rd     !<
    961        REAL(wp)     ::  rl     !<
    962        REAL(wp)     ::  rn     !<
    963        REAL(wp)     ::  rr     !<
    964        REAL(wp)     ::  rs     !<
    965        REAL(wp)     ::  rt     !<
    966946       REAL(wp)     ::  u_comp !<
    967947       REAL(wp)     ::  v_comp !<
     
    12321212                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    12331213                                         )
    1234 !
    1235 !--       Apply monotonic adjustment.
    1236           IF ( monotonic_adjustment )  THEN
    1237 !
    1238 !--          At first, calculate first order fluxes.
    1239              u_comp = u(k,j,i) - u_gtrans
    1240              fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )               &
    1241                    -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )               &
    1242                        ) * adv_sca_1
    1243 
    1244              u_comp = u(k,j,i+1) - u_gtrans
    1245              fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )               &
    1246                    -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )               &
    1247                        ) * adv_sca_1
    1248 
    1249              v_comp = v(k,j,i) - v_gtrans
    1250              fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )               &
    1251                    -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )               &
    1252                        ) * adv_sca_1
    1253 
    1254              v_comp = v(k,j+1,i) - v_gtrans
    1255              fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )               &
    1256                    -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )               &
    1257                        ) * adv_sca_1
    1258 
    1259              fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
    1260                    -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
    1261                        ) * adv_sca_1 * rho_air_zw(k)
    1262 
    1263              ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
    1264                    -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
    1265                       ) * adv_sca_1 * rho_air_zw(k)
    1266 !
    1267 !--          Calculate ratio of upwind gradients. Note, Min/Max is just to
    1268 !--          avoid if statements.
    1269              rl     = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) *                  &
    1270                            ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /    &
    1271                                 ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )      &
    1272                               ) +                                             &
    1273                         MIN( 0.0_wp, u(k,j,i) - u_gtrans ) *                  &
    1274                            ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /    &
    1275                                 ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )      &
    1276                               )                                               &
    1277                       ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
    1278 
    1279              rr     = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) *                &
    1280                            ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /    &
    1281                                 ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )      &
    1282                               ) +                                             &
    1283                         MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) *                &
    1284                            ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /    &
    1285                                 ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )      &
    1286                               )                                               &
    1287                       ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
    1288 
    1289              rs     = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) *                  &
    1290                            ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /    &
    1291                                 ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )      &
    1292                               ) +                                             &
    1293                         MIN( 0.0_wp, v(k,j,i) - v_gtrans ) *                  &
    1294                            ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /    &
    1295                                 ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )      &
    1296                               )                                               &
    1297                       ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
    1298 
    1299              rn     = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) *                &
    1300                            ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /    &
    1301                                 ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )      &
    1302                               ) +                                             &
    1303                        MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) *                 &
    1304                            ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /    &
    1305                                 ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )      &
    1306                               )                                               &
    1307                       ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
    1308 !   
    1309 !--          Reuse k_mm and compute k_mmm for the vertical gradient ratios.
    1310 !--          Note, for vertical advection below the third grid point above
    1311 !--          surface ( or below the model top) rd and rt are set to 0, i.e.
    1312 !--          use of first order scheme is enforced.
    1313              k_mmm  = k - 3 * ibit8
    1314 
    1315              rd     = ( MAX( 0.0_wp, w(k-1,j,i) ) *                           &
    1316                            ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) /  &
    1317                                 ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )    &
    1318                               ) +                                             &
    1319                         MIN( 0.0_wp, w(k-1,j,i) ) *                           &
    1320                            ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /      &
    1321                                 ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )   &
    1322                               )                                               &
    1323                       ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
    1324  
    1325              rt     = ( MAX( 0.0_wp, w(k,j,i) ) *                             &
    1326                            ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /    &
    1327                                 ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )      &
    1328                               ) +                                             &
    1329                         MIN( 0.0_wp, w(k,j,i) ) *                             &
    1330                            ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /    &
    1331                                 ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )      &
    1332                               )                                               &
    1333                       ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
    1334 !
    1335 !--           Calculate empirical limiter function (van Albada2 limiter).
    1336               phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) /                   &
    1337                                          ( rl**2 + 1.0_wp ) )
    1338               phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) /                   &
    1339                                          ( rr**2 + 1.0_wp ) )
    1340               phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) /                   &
    1341                                          ( rs**2 + 1.0_wp ) )
    1342               phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) /                   &
    1343                                          ( rn**2 + 1.0_wp ) )
    1344               phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) /                   &
    1345                                          ( rd**2 + 1.0_wp ) )
    1346               phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) /                   &
    1347                                          ( rt**2 + 1.0_wp ) )
    1348 !
    1349 !--           Calculate the resulting monotone flux.
    1350               swap_flux_x_local(k,j,tn) = fl_1 - phi_l *                      &
    1351                                         ( fl_1 - swap_flux_x_local(k,j,tn)  )
    1352               flux_r(k)                 = fr_1 - phi_r *                      &
    1353                                         ( fr_1 - flux_r(k)                  )
    1354               swap_flux_y_local(k,tn)   = fs_1 - phi_s *                      &
    1355                                         ( fs_1 - swap_flux_y_local(k,tn)    )
    1356               flux_n(k)                 = fn_1 - phi_n *                      &
    1357                                         ( fn_1 - flux_n(k)                  )
    1358               flux_d                    = fd_1 - phi_d *                      &
    1359                                         ( fd_1 - flux_d                     )
    1360               flux_t(k)                 = ft_1 - phi_t *                      &
    1361                                         ( ft_1 - flux_t(k)                  )
    1362 !
    1363 !--          Moreover, modify dissipation flux according to the limiter.
    1364              swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l
    1365              diss_r(k)                 = diss_r(k)                 * phi_r
    1366              swap_diss_y_local(k,tn)   = swap_diss_y_local(k,tn)   * phi_s
    1367              diss_n(k)                 = diss_n(k)                 * phi_n
    1368              diss_d                    = diss_d                    * phi_d
    1369              diss_t(k)                 = diss_t(k)                 * phi_t
    1370 
    1371           ENDIF
    13721214!
    13731215!--       Calculate the divergence of the velocity field. A respective
     
    14781320                             ( sk(k_ppp,j,i) - sk(k_mm,j,i) )                 &
    14791321                                         )
    1480 
    1481 
    1482 !
    1483 !--       Apply monotonic adjustment.
    1484           IF ( monotonic_adjustment )  THEN
    1485 !
    1486 !--          At first, calculate first order fluxes.
    1487              u_comp = u(k,j,i) - u_gtrans
    1488              fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )               &
    1489                    -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )               &
    1490                        ) * adv_sca_1
    1491 
    1492              u_comp = u(k,j,i+1) - u_gtrans
    1493              fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )               &
    1494                    -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )               &
    1495                        ) * adv_sca_1
    1496 
    1497              v_comp = v(k,j,i) - v_gtrans
    1498              fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )               &
    1499                    -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )               &
    1500                        ) * adv_sca_1
    1501 
    1502              v_comp = v(k,j+1,i) - v_gtrans
    1503              fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )               &
    1504                    -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )               &
    1505                        ) * adv_sca_1
    1506 
    1507              fd_1   =  ( w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )           &
    1508                    -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )           &
    1509                        ) * adv_sca_1  * rho_air_zw(k)
    1510 
    1511              ft_1   =  ( w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )             &
    1512                    -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )             &
    1513                        ) * adv_sca_1  * rho_air_zw(k)
    1514 !
    1515 !--          Calculate ratio of upwind gradients. Note, Min/Max is just to
    1516 !--          avoid if statements.
    1517              rl     = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) *                  &
    1518                            ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /    &
    1519                                 ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )      &
    1520                               ) +                                             &
    1521                         MIN( 0.0_wp, u(k,j,i) - u_gtrans ) *                  &
    1522                            ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /    &
    1523                                 ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )      &
    1524                               )                                               &
    1525                       ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
    1526 
    1527              rr     = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) *                &
    1528                            ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /    &
    1529                                 ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )      &
    1530                               ) +                                             &
    1531                         MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) *                &
    1532                            ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /    &
    1533                                 ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )      &
    1534                               )                                               &
    1535                       ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
    1536 
    1537              rs     = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) *                  &
    1538                            ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /    &
    1539                                 ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )      &
    1540                               ) +                                             &
    1541                         MIN( 0.0_wp, v(k,j,i) - v_gtrans ) *                  &
    1542                            ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /    &
    1543                                 ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )      &
    1544                               )                                               &
    1545                       ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
    1546 
    1547              rn     = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) *                &
    1548                            ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /    &
    1549                                 ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )      &
    1550                               ) +                                             &
    1551                        MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) *                 &
    1552                            ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /    &
    1553                                 ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )      &
    1554                               )                                               &
    1555                       ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
    1556 !
    1557 !--          Reuse k_mm and compute k_mmm for the vertical gradient ratios.
    1558 !--          Note, for vertical advection below the third grid point above
    1559 !--          surface ( or below the model top) rd and rt are set to 0, i.e.
    1560 !--          use of first order scheme is enforced.
    1561              k_mmm  = k - 3 * ibit8
    1562 
    1563              rd     = ( MAX( 0.0_wp, w(k-1,j,i) ) *                           &
    1564                            ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) /  &
    1565                                 ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )    &
    1566                               ) +                                             &
    1567                         MIN( 0.0_wp, w(k-1,j,i) ) *                           &
    1568                            ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /      &
    1569                                 ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )   &
    1570                               )                                               &
    1571                       ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
    1572  
    1573              rt     = ( MAX( 0.0_wp, w(k,j,i) ) *                             &
    1574                            ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /    &
    1575                                 ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )      &
    1576                               ) +                                             &
    1577                         MIN( 0.0_wp, w(k,j,i) ) *                             &
    1578                            ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /    &
    1579                                 ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )      &
    1580                               )                                               &
    1581                       ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
    1582 !
    1583 !--           Calculate empirical limiter function (van Albada2 limiter).
    1584               phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) /                   &
    1585                                          ( rl**2 + 1.0_wp ) )
    1586               phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) /                   &
    1587                                          ( rr**2 + 1.0_wp ) )
    1588               phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) /                   &
    1589                                          ( rs**2 + 1.0_wp ) )
    1590               phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) /                   &
    1591                                          ( rn**2 + 1.0_wp ) )
    1592               phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) /                   &
    1593                                          ( rd**2 + 1.0_wp ) )
    1594               phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) /                   &
    1595                                          ( rt**2 + 1.0_wp ) )
    1596 !
    1597 !--           Calculate the resulting monotone flux.
    1598               swap_flux_x_local(k,j,tn) = fl_1 - phi_l *                      &
    1599                                         ( fl_1 - swap_flux_x_local(k,j,tn)  )
    1600               flux_r(k)                 = fr_1 - phi_r *                      &
    1601                                         ( fr_1 - flux_r(k)                  )
    1602               swap_flux_y_local(k,tn)   = fs_1 - phi_s *                      &
    1603                                         ( fs_1 - swap_flux_y_local(k,tn)    )
    1604               flux_n(k)                 = fn_1 - phi_n *                      &
    1605                                         ( fn_1 - flux_n(k)                  )
    1606               flux_d                    = fd_1 - phi_d *                      &
    1607                                         ( fd_1 - flux_d                     )
    1608               flux_t(k)                 = ft_1 - phi_t *                      &
    1609                                         ( ft_1 - flux_t(k)                  )
    1610 !
    1611 !--          Moreover, modify dissipation flux according to the limiter.
    1612              swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l
    1613              diss_r(k)                 = diss_r(k)                 * phi_r
    1614              swap_diss_y_local(k,tn)   = swap_diss_y_local(k,tn)   * phi_s
    1615              diss_n(k)                 = diss_n(k)                 * phi_n
    1616              diss_d                    = diss_d                    * phi_d
    1617              diss_t(k)                 = diss_t(k)                 * phi_t
    1618 
    1619           ENDIF
    16201322!
    16211323!--       Calculate the divergence of the velocity field. A respective
     
    32072909
    32082910       USE control_parameters,                                                &
    3209            ONLY:  intermediate_timestep_count, monotonic_adjustment, u_gtrans,&
    3210                   v_gtrans
     2911           ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans
    32112912
    32122913       USE grid_variables,                                                    &
     
    32552956       REAL(wp) ::  div    !<
    32562957       REAL(wp) ::  flux_d !<
    3257        REAL(wp) ::  fd_1   !<
    3258        REAL(wp) ::  fl_1   !<
    3259        REAL(wp) ::  fn_1   !<
    3260        REAL(wp) ::  fr_1   !<
    3261        REAL(wp) ::  fs_1   !<
    3262        REAL(wp) ::  ft_1   !<
    3263        REAL(wp) ::  phi_d  !<
    3264        REAL(wp) ::  phi_l  !<
    3265        REAL(wp) ::  phi_n  !<
    3266        REAL(wp) ::  phi_r  !<
    3267        REAL(wp) ::  phi_s  !<
    3268        REAL(wp) ::  phi_t  !<
    3269        REAL(wp) ::  rd     !<
    3270        REAL(wp) ::  rl     !<
    3271        REAL(wp) ::  rn     !<
    3272        REAL(wp) ::  rr     !<
    3273        REAL(wp) ::  rs     !<
    3274        REAL(wp) ::  rt     !<
    32752958       REAL(wp) ::  u_comp !<
    32762959       REAL(wp) ::  v_comp !<
     
    35303213                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
    35313214                                               )
    3532 !
    3533 !--             Apply monotonic adjustment.
    3534                 IF ( monotonic_adjustment )  THEN
    3535 !
    3536 !--                At first, calculate first order fluxes.
    3537                    u_comp = u(k,j,i) - u_gtrans
    3538                    fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
    3539                          -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
    3540                              ) * adv_sca_1
    3541 
    3542                    u_comp = u(k,j,i+1) - u_gtrans
    3543                    fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
    3544                          -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
    3545                              ) * adv_sca_1
    3546 
    3547                    v_comp = v(k,j,i) - v_gtrans
    3548                    fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
    3549                          -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
    3550                              ) * adv_sca_1
    3551 
    3552                    v_comp = v(k,j+1,i) - v_gtrans
    3553                    fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
    3554                          -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
    3555                              ) * adv_sca_1
    3556 
    3557                    fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
    3558                          -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
    3559                             ) * adv_sca_1 * rho_air_zw(k)
    3560 
    3561                    ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
    3562                          -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
    3563                             ) * adv_sca_1 * rho_air_zw(k)
    3564 !
    3565 !--                Calculate ratio of upwind gradients. Note, Min/Max is just
    3566 !--                to avoid if statements.
    3567                    rl     = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) *            &
    3568                                ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
    3569                                     ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
    3570                                   ) +                                         &
    3571                               MIN( 0.0_wp, u(k,j,i) - u_gtrans ) *            &
    3572                                ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
    3573                                     ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
    3574                                   )                                           &
    3575                             ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
    3576 
    3577                    rr     = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) *          &
    3578                                ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
    3579                                     ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
    3580                                   ) +                                         &
    3581                               MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) *          &
    3582                                ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
    3583                                     ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
    3584                                   )                                           &
    3585                             ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
    3586 
    3587                    rs     = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) *            &
    3588                                ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
    3589                                     ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
    3590                                   ) +                                         &
    3591                               MIN( 0.0_wp, v(k,j,i) - v_gtrans ) *            &
    3592                                ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
    3593                                     ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
    3594                                   )                                           &
    3595                             ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
    3596 
    3597                    rn     = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) *          &
    3598                                ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
    3599                                     ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
    3600                                   ) +                                         &
    3601                               MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) *          &
    3602                                ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
    3603                                     ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
    3604                                   )                                           &
    3605                             ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
    3606 !   
    3607 !--                Reuse k_mm and compute k_mmm for the vertical gradient ratios.
    3608 !--                Note, for vertical advection below the third grid point above
    3609 !--                surface ( or below the model top) rd and rt are set to 0, i.e.
    3610 !--                use of first order scheme is enforced.
    3611                    k_mmm  = k - 3 * ibit8
    3612 
    3613                    rd     = ( MAX( 0.0_wp, w(k-1,j,i) ) *                     &
    3614                             ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
    3615                                  ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
    3616                                ) +                                            &
    3617                               MIN( 0.0_wp, w(k-1,j,i) ) *                     &
    3618                             ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
    3619                                  ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
    3620                                )                                              &
    3621                             ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
    3622  
    3623                    rt     = ( MAX( 0.0_wp, w(k,j,i) ) *                       &
    3624                             ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
    3625                                  ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
    3626                                ) +                                            &
    3627                               MIN( 0.0_wp, w(k,j,i) ) *                       &
    3628                             ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
    3629                                  ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
    3630                                )                                              &
    3631                             ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
    3632 !
    3633 !--                Calculate empirical limiter function (van Albada2 limiter).
    3634                    phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) /              &
    3635                                         ( rl**2 + 1.0_wp ) )
    3636                    phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) /              &
    3637                                         ( rr**2 + 1.0_wp ) )
    3638                    phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) /              &
    3639                                         ( rs**2 + 1.0_wp ) )
    3640                    phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) /              &
    3641                                         ( rn**2 + 1.0_wp ) )
    3642                    phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) /              &
    3643                                         ( rd**2 + 1.0_wp ) )
    3644                    phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) /              &
    3645                                         ( rt**2 + 1.0_wp ) )
    3646 !
    3647 !--                Calculate the resulting monotone flux.
    3648                    swap_flux_x_local(k,j) = fl_1 - phi_l *                    &
    3649                                           ( fl_1 - swap_flux_x_local(k,j) )
    3650                    flux_r(k)              = fr_1 - phi_r *                    &
    3651                                           ( fr_1 - flux_r(k)              )
    3652                    swap_flux_y_local(k)   = fs_1 - phi_s *                    &
    3653                                           ( fs_1 - swap_flux_y_local(k)   )
    3654                    flux_n(k)              = fn_1 - phi_n *                    &
    3655                                           ( fn_1 - flux_n(k)              )
    3656                    flux_d                 = fd_1 - phi_d *                    &
    3657                                           ( fd_1 - flux_d                 )
    3658                    flux_t(k)              = ft_1 - phi_t *                    &
    3659                                           ( ft_1 - flux_t(k)              )
    3660 !
    3661 !--                Moreover, modify dissipation flux according to the limiter.
    3662                    swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l
    3663                    diss_r(k)              = diss_r(k)              * phi_r
    3664                    swap_diss_y_local(k)   = swap_diss_y_local(k)   * phi_s
    3665                    diss_n(k)              = diss_n(k)              * phi_n
    3666                    diss_d                 = diss_d                 * phi_d
    3667                    diss_t(k)              = diss_t(k)              * phi_t
    3668 
    3669                 ENDIF
    36703215!
    36713216!--             Calculate the divergence of the velocity field. A respective
     
    37743319                                   ( sk(k_ppp,j,i) - sk(k_mm,j,i) )           &
    37753320                                               )
    3776 !
    3777 !--             Apply monotonic adjustment.
    3778                 IF ( monotonic_adjustment )  THEN
    3779 !
    3780 !--                At first, calculate first order fluxes.
    3781                    u_comp = u(k,j,i) - u_gtrans
    3782                    fl_1   =  ( u_comp   * ( sk(k,j,i) + sk(k,j,i-1) )         &
    3783                          -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) )         &
    3784                              ) * adv_sca_1
    3785 
    3786                    u_comp = u(k,j,i+1) - u_gtrans
    3787                    fr_1   =  ( u_comp   * ( sk(k,j,i+1) + sk(k,j,i) )         &
    3788                          -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) )         &
    3789                              ) * adv_sca_1
    3790 
    3791                    v_comp = v(k,j,i) - v_gtrans
    3792                    fs_1   =  ( v_comp   * ( sk(k,j,i) + sk(k,j-1,i) )         &
    3793                          -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) )         &
    3794                              ) * adv_sca_1
    3795 
    3796                    v_comp = v(k,j+1,i) - v_gtrans
    3797                    fn_1   =  ( v_comp   * ( sk(k,j+1,i) + sk(k,j,i) )         &
    3798                          -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) )         &
    3799                              ) * adv_sca_1
    3800 
    3801                    fd_1   = (  w(k-1,j,i)   * ( sk(k,j,i) + sk(k-1,j,i) )     &
    3802                          -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) )     &
    3803                             ) * adv_sca_1 * rho_air_zw(k)
    3804 
    3805                    ft_1   = (  w(k,j,i)   * ( sk(k+1,j,i) + sk(k,j,i) )       &
    3806                          -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) )       &
    3807                             ) * adv_sca_1 * rho_air_zw(k)
    3808 !
    3809 !--                Calculate ratio of upwind gradients. Note, Min/Max is just
    3810 !--                to avoid if statements.
    3811                    rl     = ( MAX( 0.0_wp, u(k,j,i) - u_gtrans ) *            &
    3812                                ABS( ( sk(k,j,i-1) - sk(k,j,i-2)            ) /&
    3813                                     ( sk(k,j,i)   - sk(k,j,i-1) + 1E-20_wp )  &
    3814                                   ) +                                         &
    3815                               MIN( 0.0_wp, u(k,j,i) - u_gtrans ) *            &
    3816                                ABS( ( sk(k,j,i)   - sk(k,j,i+1)            ) /&
    3817                                     ( sk(k,j,i-1) - sk(k,j,i)   + 1E-20_wp )  &
    3818                                   )                                           &
    3819                             ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp )
    3820 
    3821                    rr     = ( MAX( 0.0_wp, u(k,j,i+1) - u_gtrans ) *          &
    3822                                ABS( ( sk(k,j,i)   - sk(k,j,i-1)            ) /&
    3823                                     ( sk(k,j,i+1) - sk(k,j,i)   + 1E-20_wp )  &
    3824                                   ) +                                         &
    3825                               MIN( 0.0_wp, u(k,j,i+1) - u_gtrans ) *          &
    3826                                ABS( ( sk(k,j,i+1) - sk(k,j,i+2)            ) /&
    3827                                     ( sk(k,j,i)   - sk(k,j,i+1) + 1E-20_wp )  &
    3828                                   )                                           &
    3829                             ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp )
    3830 
    3831                    rs     = ( MAX( 0.0_wp, v(k,j,i) - v_gtrans ) *            &
    3832                                ABS( ( sk(k,j-1,i) - sk(k,j-2,i)            ) /&
    3833                                     ( sk(k,j,i)   - sk(k,j-1,i) + 1E-20_wp )  &
    3834                                   ) +                                         &
    3835                               MIN( 0.0_wp, v(k,j,i) - v_gtrans ) *            &
    3836                                ABS( ( sk(k,j,i)   - sk(k,j+1,i)            ) /&
    3837                                     ( sk(k,j-1,i) - sk(k,j,i)   + 1E-20_wp )  &
    3838                                   )                                           &
    3839                             ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp )
    3840 
    3841                    rn     = ( MAX( 0.0_wp, v(k,j+1,i) - v_gtrans ) *          &
    3842                                ABS( ( sk(k,j,i)   - sk(k,j-1,i)            ) /&
    3843                                     ( sk(k,j+1,i) - sk(k,j,i)   + 1E-20_wp )  &
    3844                                   ) +                                         &
    3845                               MIN( 0.0_wp, v(k,j+1,i) - v_gtrans ) *          &
    3846                                ABS( ( sk(k,j+1,i) - sk(k,j+2,i)            ) /&
    3847                                     ( sk(k,j,i)   - sk(k,j+1,i) + 1E-20_wp )  &
    3848                                   )                                           &
    3849                             ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp )     
    3850 !   
    3851 !--                Reuse k_mm and compute k_mmm for the vertical gradient ratios.
    3852 !--                Note, for vertical advection below the third grid point above
    3853 !--                surface ( or below the model top) rd and rt are set to 0, i.e.
    3854 !--                use of first order scheme is enforced.
    3855                    k_mmm  = k - 3 * ibit8
    3856 
    3857                    rd     = ( MAX( 0.0_wp, w(k-1,j,i) ) *                     &
    3858                             ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i)           ) / &
    3859                                  ( sk(k-1,j,i)  - sk(k_mm,j,i) + 1E-20_wp )   &
    3860                                ) +                                            &
    3861                               MIN( 0.0_wp, w(k-1,j,i) ) *                     &
    3862                             ABS( ( sk(k-1,j,i) - sk(k,j,i)            ) /     &
    3863                                  ( sk(k_mm,j,i) - sk(k-1,j,i)   + 1E-20_wp )  &
    3864                                )                                              &
    3865                             ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp )
    3866  
    3867                    rt     = ( MAX( 0.0_wp, w(k,j,i) ) *                       &
    3868                             ABS( ( sk(k,j,i)   - sk(k-1,j,i)            ) /   &
    3869                                  ( sk(k+1,j,i) - sk(k,j,i)   + 1E-20_wp )     &
    3870                                ) +                                            &
    3871                               MIN( 0.0_wp, w(k,j,i) ) *                       &
    3872                             ABS( ( sk(k+1,j,i) - sk(k_pp,j,i)           ) /   &
    3873                                  ( sk(k,j,i)   - sk(k+1,j,i) + 1E-20_wp )     &
    3874                                )                                              &
    3875                             ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp )
    3876 !
    3877 !--                Calculate empirical limiter function (van Albada2 limiter).
    3878                    phi_l = MIN( 1.0_wp, ( 2.0_wp * ABS( rl ) ) /              &
    3879                                         ( rl**2 + 1.0_wp ) )
    3880                    phi_r = MIN( 1.0_wp, ( 2.0_wp * ABS( rr ) ) /              &
    3881                                         ( rr**2 + 1.0_wp ) )
    3882                    phi_s = MIN( 1.0_wp, ( 2.0_wp * ABS( rs ) ) /              &
    3883                                         ( rs**2 + 1.0_wp ) )
    3884                    phi_n = MIN( 1.0_wp, ( 2.0_wp * ABS( rn ) ) /              &
    3885                                         ( rn**2 + 1.0_wp ) )
    3886                    phi_d = MIN( 1.0_wp, ( 2.0_wp * ABS( rd ) ) /              &
    3887                                         ( rd**2 + 1.0_wp ) )
    3888                    phi_t = MIN( 1.0_wp, ( 2.0_wp * ABS( rt ) ) /              &
    3889                                         ( rt**2 + 1.0_wp ) )
    3890 !
    3891 !--                Calculate the resulting monotone flux.
    3892                    swap_flux_x_local(k,j) = fl_1 - phi_l *                    &
    3893                                           ( fl_1 - swap_flux_x_local(k,j) )
    3894                    flux_r(k)              = fr_1 - phi_r *                    &
    3895                                           ( fr_1 - flux_r(k)              )
    3896                    swap_flux_y_local(k)   = fs_1 - phi_s *                    &
    3897                                           ( fs_1 - swap_flux_y_local(k)   )
    3898                    flux_n(k)              = fn_1 - phi_n *                    &
    3899                                           ( fn_1 - flux_n(k)              )
    3900                    flux_d                 = fd_1 - phi_d *                    &
    3901                                           ( fd_1 - flux_d                 )
    3902                    flux_t(k)              = ft_1 - phi_t *                    &
    3903                                           ( ft_1 - flux_t(k)              )
    3904 !
    3905 !--                Moreover, modify dissipation flux according to the limiter.
    3906                    swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l
    3907                    diss_r(k)              = diss_r(k)              * phi_r
    3908                    swap_diss_y_local(k)   = swap_diss_y_local(k)   * phi_s
    3909                    diss_n(k)              = diss_n(k)              * phi_n
    3910                    diss_d                 = diss_d                 * phi_d
    3911                    diss_t(k)              = diss_t(k)              * phi_t
    3912 
    3913                 ENDIF
    39143321!
    39153322!--             Calculate the divergence of the velocity field. A respective
Note: See TracChangeset for help on using the changeset viewer.