Ignore:
Timestamp:
Jul 12, 2016 4:34:24 PM (5 years ago)
Author:
suehring
Message:

Separate balance equations for humidity and passive_scalar

File:
1 edited

Legend:

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

    r1917 r1960  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    192192    USE arrays_3d,                                                             &
    193193        ONLY:  diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr,            &
    194                diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,            &
    195                diss_s_qr, diss_s_sa, e, e_p, flux_s_e, flux_s_nr, flux_s_pt,   &
    196                flux_s_q, flux_s_qr, flux_s_sa, flux_l_e, flux_l_nr,            &
    197                flux_l_pt, flux_l_q, flux_l_qr, flux_l_sa, nr, nr_p, nrsws,     &
    198                nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p, prho, q, q_init,     &
    199                q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf, rdf_sc,         &
    200                ref_state, rho, sa, sa_init, sa_p, saswsb, saswst, shf, tend,   &
    201                te_m, tnr_m, tpt_m, tq_m, tqr_m, tsa_m, tswst, tu_m, tv_m,      &
     194               diss_l_s, diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,  &
     195               diss_s_qr, diss_s_s, diss_s_sa, e, e_p, flux_s_e, flux_s_nr,    &
     196               flux_s_pt, flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e,  &
     197               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s, flux_l_sa, &
     198               nr, nr_p, nrsws, nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p,     &
     199               prho, q, q_init, q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf,&
     200               rdf_sc, ref_state, rho, s, s_init, s_p, sa, sa_init, sa_p,      &
     201               saswsb, saswst, shf, ssws, sswst, tend,                         &
     202               te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tswst, tu_m, tv_m,&
    202203               tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
    203204       
     
    216217               use_upstream_for_tke, wall_heatflux,                            &
    217218               wall_nrflux, wall_qflux, wall_qflux, wall_qflux, wall_qrflux,   &
    218                wall_salinityflux, ws_scheme_mom, ws_scheme_sca
     219               wall_salinityflux, wall_sflux, ws_scheme_mom, ws_scheme_sca
    219220
    220221    USE cpulog,                                                                &
     
    724725
    725726!
    726 !--       If required, compute prognostic equation for total water content /
    727 !--       scalar
    728           IF ( humidity  .OR.  passive_scalar )  THEN
     727!--       If required, compute prognostic equation for total water content
     728          IF ( humidity )  THEN
    729729
    730730!
     
    745745
    746746!
    747 !--          Sink or source of scalar concentration due to canopy elements
     747!--          Sink or source of humidity due to canopy elements
    748748             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
    749749
     
    881881
    882882          ENDIF
    883 
     883         
     884!
     885!--       If required, compute prognostic equation for scalar
     886          IF ( passive_scalar )  THEN
     887!
     888!--          Tendency-terms for total water content / scalar
     889             tend(:,j,i) = 0.0_wp
     890             IF ( timestep_scheme(1:5) == 'runge' ) &
     891             THEN
     892                IF ( ws_scheme_sca )  THEN
     893                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
     894                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
     895                ELSE
     896                   CALL advec_s_pw( i, j, s )
     897                ENDIF
     898             ELSE
     899                CALL advec_s_up( i, j, s )
     900             ENDIF
     901             CALL diffusion_s( i, j, s, ssws, sswst, wall_sflux )
     902
     903!
     904!--          Sink or source of scalar concentration due to canopy elements
     905             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
     906
     907!
     908!--          Large scale advection, still need to be extended for scalars
     909!              IF ( large_scale_forcing )  THEN
     910!                 CALL ls_advec( i, j, simulated_time, 's' )
     911!              ENDIF
     912
     913!
     914!--          Nudging, still need to be extended for scalars
     915!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
     916
     917!
     918!--          If required compute influence of large-scale subsidence/ascent.
     919!--          Note, the last argument is of no meaning in this case, as it is
     920!--          only used in conjunction with large_scale_forcing, which is to
     921!--          date not implemented for scalars.
     922             IF ( large_scale_subsidence  .AND.                                &
     923                  .NOT. use_subsidence_tendencies  .AND.                       &
     924                  .NOT. large_scale_forcing )  THEN
     925                CALL subsidence( i, j, tend, s, s_init, 3 )
     926             ENDIF
     927
     928             CALL user_actions( i, j, 's-tendency' )
     929
     930!
     931!--          Prognostic equation for scalar
     932             DO  k = nzb_s_inner(j,i)+1, nzt
     933                s_p(k,j,i) = s(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
     934                                                  tsc(3) * ts_m(k,j,i) )       &
     935                                      - tsc(5) * rdf_sc(k) *                   &
     936                                        ( s(k,j,i) - s_init(k) )
     937                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
     938             ENDDO
     939
     940!
     941!--          Calculate tendencies for the next Runge-Kutta step
     942             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     943                IF ( intermediate_timestep_count == 1 )  THEN
     944                   DO  k = nzb_s_inner(j,i)+1, nzt
     945                      ts_m(k,j,i) = tend(k,j,i)
     946                   ENDDO
     947                ELSEIF ( intermediate_timestep_count < &
     948                         intermediate_timestep_count_max )  THEN
     949                   DO  k = nzb_s_inner(j,i)+1, nzt
     950                      ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
     951                                     5.3125_wp * ts_m(k,j,i)
     952                   ENDDO
     953                ENDIF
     954             ENDIF
     955
     956          ENDIF         
    884957!
    885958!--       If required, compute prognostic equation for turbulent kinetic
     
    14321505
    14331506!
    1434 !-- If required, compute prognostic equation for total water content / scalar
    1435     IF ( humidity  .OR.  passive_scalar )  THEN
    1436 
    1437        CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
     1507!-- If required, compute prognostic equation for total water content
     1508    IF ( humidity )  THEN
     1509
     1510       CALL cpu_log( log_point(29), 'q-equation', 'start' )
    14381511
    14391512!
     
    14701543       
    14711544!
    1472 !--    Sink or source of scalar concentration due to canopy elements
     1545!--    Sink or source of humidity due to canopy elements
    14731546       IF ( plant_canopy ) CALL pcm_tendency( 5 )
    14741547
     
    14931566
    14941567!
    1495 !--    Prognostic equation for total water content / scalar
     1568!--    Prognostic equation for total water content
    14961569       DO  i = nxl, nxr
    14971570          DO  j = nys, nyn
     
    15291602       ENDIF
    15301603
    1531        CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
     1604       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
    15321605
    15331606!
     
    16861759
    16871760    ENDIF
    1688 
     1761!
     1762!-- If required, compute prognostic equation for scalar
     1763    IF ( passive_scalar )  THEN
     1764
     1765       CALL cpu_log( log_point(66), 's-equation', 'start' )
     1766
     1767!
     1768!--    Scalar/q-tendency terms with communication
     1769       sbt = tsc(2)
     1770       IF ( scalar_advec == 'bc-scheme' )  THEN
     1771
     1772          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1773!
     1774!--          Bott-Chlond scheme always uses Euler time step. Thus:
     1775             sbt = 1.0_wp
     1776          ENDIF
     1777          tend = 0.0_wp
     1778          CALL advec_s_bc( s, 's' )
     1779
     1780       ENDIF
     1781
     1782!
     1783!--    Scalar/q-tendency terms with no communication
     1784       IF ( scalar_advec /= 'bc-scheme' )  THEN
     1785          tend = 0.0_wp
     1786          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1787             IF ( ws_scheme_sca )  THEN
     1788                CALL advec_s_ws( s, 's' )
     1789             ELSE
     1790                CALL advec_s_pw( s )
     1791             ENDIF
     1792          ELSE
     1793             CALL advec_s_up( s )
     1794          ENDIF
     1795       ENDIF
     1796
     1797       CALL diffusion_s( s, ssws, sswst, wall_sflux )
     1798       
     1799!
     1800!--    Sink or source of humidity due to canopy elements
     1801       IF ( plant_canopy ) CALL pcm_tendency( 7 )
     1802
     1803!
     1804!--    Large scale advection. Not implemented for scalars so far.
     1805!        IF ( large_scale_forcing )  THEN
     1806!           CALL ls_advec( simulated_time, 'q' )
     1807!        ENDIF
     1808
     1809!
     1810!--    Nudging. Not implemented for scalars so far.
     1811!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
     1812
     1813!
     1814!--    If required compute influence of large-scale subsidence/ascent.
     1815!--    Not implemented for scalars so far.
     1816       IF ( large_scale_subsidence  .AND.                                      &
     1817            .NOT. use_subsidence_tendencies  .AND.                             &
     1818            .NOT. large_scale_forcing )  THEN
     1819         CALL subsidence( tend, s, s_init, 3 )
     1820       ENDIF
     1821
     1822       CALL user_actions( 's-tendency' )
     1823
     1824!
     1825!--    Prognostic equation for total water content
     1826       DO  i = nxl, nxr
     1827          DO  j = nys, nyn
     1828             DO  k = nzb_s_inner(j,i)+1, nzt
     1829                s_p(k,j,i) = s(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
     1830                                                  tsc(3) * ts_m(k,j,i) )       &
     1831                                      - tsc(5) * rdf_sc(k) *                   &
     1832                                        ( s(k,j,i) - s_init(k) )
     1833                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
     1834             ENDDO
     1835          ENDDO
     1836       ENDDO
     1837
     1838!
     1839!--    Calculate tendencies for the next Runge-Kutta step
     1840       IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1841          IF ( intermediate_timestep_count == 1 )  THEN
     1842             DO  i = nxl, nxr
     1843                DO  j = nys, nyn
     1844                   DO  k = nzb_s_inner(j,i)+1, nzt
     1845                      ts_m(k,j,i) = tend(k,j,i)
     1846                   ENDDO
     1847                ENDDO
     1848             ENDDO
     1849          ELSEIF ( intermediate_timestep_count < &
     1850                   intermediate_timestep_count_max )  THEN
     1851             DO  i = nxl, nxr
     1852                DO  j = nys, nyn
     1853                   DO  k = nzb_s_inner(j,i)+1, nzt
     1854                      ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * ts_m(k,j,i)
     1855                   ENDDO
     1856                ENDDO
     1857             ENDDO
     1858          ENDIF
     1859       ENDIF
     1860
     1861       CALL cpu_log( log_point(66), 's-equation', 'stop' )
     1862
     1863    ENDIF
    16891864!
    16901865!-- If required, compute prognostic equation for turbulent kinetic
     
    22262401
    22272402!
    2228 !-- If required, compute prognostic equation for total water content / scalar
    2229     IF ( humidity  .OR.  passive_scalar )  THEN
    2230 
    2231        CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
     2403!-- If required, compute prognostic equation for total water content
     2404    IF ( humidity )  THEN
     2405
     2406       CALL cpu_log( log_point(29), 'q-equation', 'start' )
    22322407
    22332408!
     
    23072482       ENDDO
    23082483
    2309        CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
     2484       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
    23102485
    23112486!
     
    24302605    ENDIF
    24312606
     2607!
     2608!-- If required, compute prognostic equation for scalar
     2609    IF ( passive_scalar )  THEN
     2610
     2611       CALL cpu_log( log_point(66), 's-equation', 'start' )
     2612
     2613!
     2614!--    Scalar/q-tendency terms with communication
     2615       sbt = tsc(2)
     2616       IF ( scalar_advec == 'bc-scheme' )  THEN
     2617
     2618          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     2619!
     2620!--          Bott-Chlond scheme always uses Euler time step. Thus:
     2621             sbt = 1.0_wp
     2622          ENDIF
     2623          tend = 0.0_wp
     2624          CALL advec_s_bc( s, 's' )
     2625
     2626       ENDIF
     2627
     2628!
     2629!--    Scalar/q-tendency terms with no communication
     2630       IF ( scalar_advec /= 'bc-scheme' )  THEN
     2631          tend = 0.0_wp
     2632          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2633             IF ( ws_scheme_sca )  THEN
     2634                CALL advec_s_ws( s, 's' )
     2635             ELSE
     2636                CALL advec_s_pw( s )
     2637             ENDIF
     2638          ELSE
     2639             CALL advec_s_up( s )
     2640          ENDIF
     2641       ENDIF
     2642
     2643       CALL diffusion_s( s, ssws, sswst, wall_sflux )
     2644
     2645!
     2646!--    Sink or source of scalar concentration due to canopy elements
     2647       IF ( plant_canopy ) CALL pcm_tendency( 7 )
     2648
     2649!
     2650!--    Large scale advection. Not implemented so far.
     2651!        IF ( large_scale_forcing )  THEN
     2652!           CALL ls_advec( simulated_time, 's' )
     2653!        ENDIF
     2654
     2655!
     2656!--    Nudging. Not implemented so far.
     2657!        IF ( nudging )  CALL nudge( simulated_time, 's' )
     2658
     2659!
     2660!--    If required compute influence of large-scale subsidence/ascent.
     2661!--    Not implemented so far.
     2662       IF ( large_scale_subsidence  .AND.                                      &
     2663            .NOT. use_subsidence_tendencies  .AND.                             &
     2664            .NOT. large_scale_forcing )  THEN
     2665         CALL subsidence( tend, s, s_init, 3 )
     2666       ENDIF
     2667
     2668       CALL user_actions( 's-tendency' )
     2669
     2670!
     2671!--    Prognostic equation for total water content / scalar
     2672       DO  i = i_left, i_right
     2673          DO  j = j_south, j_north
     2674             DO  k = nzb_s_inner(j,i)+1, nzt
     2675                s_p(k,j,i) = s(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
     2676                                                  tsc(3) * ts_m(k,j,i) )       &
     2677                                      - tsc(5) * rdf_sc(k) *                   &
     2678                                        ( s(k,j,i) - s_init(k) )
     2679                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
     2680!
     2681!--             Tendencies for the next Runge-Kutta step
     2682                IF ( runge_step == 1 )  THEN
     2683                   ts_m(k,j,i) = tend(k,j,i)
     2684                ELSEIF ( runge_step == 2 )  THEN
     2685                   ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * ts_m(k,j,i)
     2686                ENDIF
     2687             ENDDO
     2688          ENDDO
     2689       ENDDO
     2690
     2691       CALL cpu_log( log_point(66), 's-equation', 'stop' )
     2692
     2693    ENDIF
    24322694!
    24332695!-- If required, compute prognostic equation for turbulent kinetic
Note: See TracChangeset for help on using the changeset viewer.