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/flow_statistics.f90

    r1919 r1960  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Separate humidity and passive scalar
    2222!
    2323! Former revisions:
     
    197197    USE arrays_3d,                                                             &
    198198        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
    199                qr, qs, qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt,  &
    200                td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us,&
    201                usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
     199               qr, qs, qsws, qswst, rho, s, sa, ss, ssws, sswst, saswsb,       &
     200               saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q,        &
     201               time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, &
     202               vswst, w, w_subs, zw
    202203       
    203204    USE cloud_parameters,                                                      &
     
    360361
    361362          DO  i = 0, threads_per_task-1
    362              sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
    363              IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
    364              IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
    365                                                    sums_wsqs_ws_l(:,i) !w*q*
     363             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i) ! w*pt*
     364             IF ( ocean          ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
     365             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)  ! w*q*
     366             IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
    366367          ENDDO
    367368
     
    440441             DO  j =  nys, nyn
    441442                DO  k = nzb_s_inner(j,i), nzt+1
    442                    sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
     443                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i) * rmask(j,i,sr)
    443444                ENDDO
    444445             ENDDO
     
    465466             ENDIF
    466467             IF ( passive_scalar )  THEN
    467                 sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
     468                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
    468469             ENDIF
    469470          ENDDO
     
    506507       IF ( passive_scalar )  THEN
    507508          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    508           CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
     509          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,       &
    509510                              MPI_REAL, MPI_SUM, comm2d, ierr )
    510511       ENDIF
     
    522523          ENDIF
    523524       ENDIF
    524        IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
     525       IF ( passive_scalar )  sums(:,117) = sums_l(:,117,0)
    525526#endif
    526527
     
    560561!
    561562!--    Passive scalar
    562        IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /                    &
    563             ngp_2dh_s_inner(:,sr)                    ! s (q)
     563       IF ( passive_scalar )  hom(:,1,117,sr) = sums(:,117) /                  &
     564            ngp_2dh_s_inner(:,sr)                    ! s
    564565
    565566!
     
    598599                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
    599600                ENDIF
    600 
     601                IF ( passive_scalar )  THEN
     602                   sums_l(k,118,tn) = sums_l(k,118,tn) + &
     603                                  ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)
     604                ENDIF
    601605!
    602606!--             Higher moments
     
    627631                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
    628632                                            qs(j,i)   * rmask(j,i,sr)
     633             ENDIF
     634             IF ( passive_scalar )  THEN
     635                sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
     636                                            ss(j,i)   * rmask(j,i,sr)
    629637             ENDIF
    630638          ENDDO
     
    738746!--             Passive scalar flux
    739747                IF ( passive_scalar )  THEN
    740                    sums_l(k,48,tn) = sums_l(k,48,tn)                           &
     748                   sums_l(k,119,tn) = sums_l(k,119,tn)                         &
    741749                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    742                                                * ( q(k+1,j,i) - q(k,j,i) )     &
     750                                               * ( s(k+1,j,i) - s(k,j,i) )     &
    743751                                               * ddzu(k+1) * rmask(j,i,sr)
    744752                ENDIF
     
    784792                ENDIF
    785793                IF ( passive_scalar )  THEN
    786                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
    787                                        qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
     794                   sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
     795                                       ssws(j,i) * rmask(j,i,sr) ! w"s"
    788796                ENDIF
    789797             ENDIF
     
    863871                ENDIF
    864872                IF ( passive_scalar )  THEN
    865                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
    866                                        qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
     873                   sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
     874                                       sswst(j,i) * rmask(j,i,sr) ! w"s"
    867875                ENDIF
    868876             ENDIF
     
    953961                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
    954962                     .OR. sr /= 0 ) )  THEN
    955                    pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +              &
    956                                     q(k+1,j,i) - hom(k+1,1,41,sr) )
    957                    sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *        &
     963                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +             &
     964                                    s(k+1,j,i) - hom(k+1,1,117,sr) )
     965                   sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
    958966                                                       rmask(j,i,sr)
    959967                ENDIF
     
    10111019                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
    10121020                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
     1021                                       rmask(j,i,sr)
     1022                  ENDIF
     1023                  IF ( passive_scalar )  THEN
     1024                     pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +            &
     1025                                      s(k+1,j,i) - hom(k+1,1,117,sr) )
     1026                     sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
    10131027                                       rmask(j,i,sr)
    10141028                  ENDIF
     
    12791293          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
    12801294          sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
     1295          sums(k,116)           = sums(k,116)           / ngp_2dh(sr)
     1296          sums(k,119)           = sums(k,119)           / ngp_2dh(sr)
    12811297          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
    12821298             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
     
    12871303             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
    12881304             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
    1289              sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
     1305             sums(k,118)           = sums(k,118)           / ngp_2dh_s_inner(k,sr)
     1306             sums(k,120:pr_palm-2) = sums(k,120:pr_palm-2) / ngp_2dh_s_inner(k,sr)
    12901307          ENDIF
    12911308       ENDDO
     
    12981315                                    ngp_2dh(sr)
    12991316       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
     1317                                    ngp_2dh(sr)
     1318       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
    13001319                                    ngp_2dh(sr)
    13011320!--    eges, e*
     
    14431462       hom(:,1,114,sr) = sums(:,114)            !: L
    14441463
     1464       IF ( passive_scalar )  THEN
     1465          hom(:,1,119,sr) = sums(:,119)     ! w"s"
     1466          hom(:,1,116,sr) = sums(:,116)     ! w*s*
     1467          hom(:,1,120,sr) = sums(:,119) + sums(:,116)    ! ws
     1468       ENDIF
     1469
    14451470       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
    14461471                                       ! u*, w'u', w'v', t* (in last profile)
     
    15891614       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
    15901615
     1616       IF ( passive_scalar )  THEN
     1617          ts_value(24,sr) = hom(nzb+13,1,119,sr)       ! w"s" ( to do ! )
     1618          ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr)   ! s*
     1619       ENDIF
     1620
    15911621!
    15921622!--    Collect land surface model timeseries
     
    16601690    USE arrays_3d,                                                             &
    16611691        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
    1662                qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
    1663                td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
    1664                uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
     1692               qsws, qswst, rho, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
     1693               td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts,      &
     1694               tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w,      &
     1695               w_subs, zw
    16651696                 
    16661697       
     
    18331864             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
    18341865             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
    1835              IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
    1836                                                    sums_wsqs_ws_l(:,i) !w*q*
     1866             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i) !w*q*
     1867             IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
    18371868          ENDDO
    18381869
     
    19381969       IF ( passive_scalar )  THEN
    19391970          !$OMP DO
    1940           !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
     1971          !$acc parallel loop gang present( s, rflags_invers, rmask, sums_l ) create( s1 )
    19411972          DO  k = nzb, nzt+1
    19421973             s1 = 0
     
    19441975             DO  i = nxl, nxr
    19451976                DO  j =  nys, nyn
    1946                    s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    1947                 ENDDO
    1948              ENDDO
    1949              sums_l(k,41,tn) = s1
     1977                   s1 = s1 + s(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     1978                ENDDO
     1979             ENDDO
     1980             sums_l(k,117,tn) = s1
    19501981          ENDDO
    19511982          !$acc end parallel loop
     
    19812012             IF ( passive_scalar )  THEN
    19822013                !$acc parallel present( sums_l )
    1983                 sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
     2014                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
    19842015                !$acc end parallel
    19852016             ENDIF
     
    20242055       IF ( passive_scalar )  THEN
    20252056          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    2026           CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
     2057          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,     &
    20272058                              MPI_REAL, MPI_SUM, comm2d, ierr )
    20282059       ENDIF
     
    20532084       IF ( passive_scalar )  THEN
    20542085          !$acc parallel present( sums, sums_l )
    2055           sums(:,41) = sums_l(:,41,0)
     2086          sums(:,117) = sums_l(:,117,0)
    20562087          !$acc end parallel
    20572088       ENDIF
     
    21022133       IF ( passive_scalar )  THEN
    21032134          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
    2104           sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
    2105           hom(:,1,41,sr) = sums(:,41)                ! s (q)
     2135          sums(:,117)     = sums(:,117) / ngp_2dh_s_inner(:,sr)
     2136          hom(:,1,117,sr) = sums(:,117)                ! s
    21062137          !$acc end parallel
    21072138       ENDIF
     
    22302261          ENDDO
    22312262          sums_l(nzb+12,pr_palm,tn) = s1
     2263          !$acc end parallel
     2264       ENDIF
     2265
     2266       IF ( passive_scalar )  THEN
     2267          !$acc parallel present( ss, rmask, sums_l ) create( s1 )
     2268          s1 = 0
     2269          !$acc loop vector collapse( 2 ) reduction( +: s1 )
     2270          DO  i = nxl, nxr
     2271             DO  j =  nys, nyn
     2272                s1 = s1 + ss(j,i) * rmask(j,i,sr)
     2273             ENDDO
     2274          ENDDO
     2275          sums_l(nzb+13,pr_palm,tn) = s1
    22322276          !$acc end parallel
    22332277       ENDIF
     
    24112455       IF ( passive_scalar )  THEN
    24122456
    2413           !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
     2457          !$acc parallel loop gang present( ddzu, kh, s, rflags_invers, rmask, sums_l ) create( s1 )
    24142458          DO  k = nzb, nzt_diff
    24152459             s1 = 0
     
    24182462                DO  j = nys, nyn
    24192463                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    2420                                     * ( q(k+1,j,i) - q(k,j,i) )                &
     2464                                    * ( s(k+1,j,i) - s(k,j,i) )                &
    24212465                                    * ddzu(k+1) * rmask(j,i,sr)                &
    24222466                                    * rflags_invers(j,i,k+1)
    24232467                ENDDO
    24242468             ENDDO
    2425              sums_l(k,48,tn) = s1
     2469             sums_l(k,119,tn) = s1
    24262470          ENDDO
    24272471          !$acc end parallel loop
     
    25322576
    25332577             !$OMP DO
    2534              !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
     2578             !$acc parallel present( ssws, rmask, sums_l ) create( s1 )
    25352579             s1 = 0
    25362580             !$acc loop vector collapse( 2 ) reduction( +: s1 )
    25372581             DO  i = nxl, nxr
    25382582                DO  j =  nys, nyn
    2539                    s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
    2540                 ENDDO
    2541              ENDDO
    2542              sums_l(nzb,48,tn) = s1
     2583                   s1 = s1 + ssws(j,i) * rmask(j,i,sr)  ! w"s"
     2584                ENDDO
     2585             ENDDO
     2586             sums_l(nzb,119,tn) = s1
    25432587             !$acc end parallel
    25442588
     
    26512695
    26522696             !$OMP DO
    2653              !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
     2697             !$acc parallel present( sswst, rmask, sums_l ) create( s1 )
    26542698             s1 = 0
    26552699             !$acc loop vector collapse( 2 ) reduction( +: s1 )
    26562700             DO  i = nxl, nxr
    26572701                DO  j =  nys, nyn
    2658                    s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    2659                 ENDDO
    2660              ENDDO
    2661              sums_l(nzt,48,tn) = s1
     2702                   s1 = s1 + sswst(j,i) * rmask(j,i,sr) ! w"s"
     2703                ENDDO
     2704             ENDDO
     2705             sums_l(nzt,119,tn) = s1
    26622706             !$acc end parallel
    26632707
     
    28902934       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
    28912935
    2892           !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
     2936          !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
    28932937          DO  k = nzb, nzt_diff
    28942938             s1 = 0
     
    28962940             DO  i = nxl, nxr
    28972941                DO  j = nys, nyn
    2898                    s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +          &
    2899                                         q(k+1,j,i) - hom(k+1,1,41,sr) )        &
     2942                   s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +          &
     2943                                        s(k+1,j,i) - hom(k+1,1,117,sr) )        &
    29002944                                    * w(k,j,i) * rmask(j,i,sr)                 &
    29012945                                    * rflags_invers(j,i,k+1)
     
    29813025                ENDDO
    29823026                sums_l(k,49,tn) = s1
     3027             ENDDO
     3028             !$acc end parallel loop
     3029
     3030          ENDIF
     3031
     3032          IF ( passive_scalar )  THEN
     3033
     3034             !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
     3035             DO  k = nzb, nzt_diff
     3036                s1 = 0
     3037                !$acc loop vector collapse( 2 ) reduction( +: s1 )
     3038                DO  i = nxl, nxr
     3039                   DO  j = nys, nyn
     3040                      s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +      &
     3041                                           s(k+1,j,i) - hom(k+1,1,117,sr) )    &
     3042                                       * w(k,j,i) * rmask(j,i,sr)              &
     3043                                       * rflags_invers(j,i,k+1)
     3044                   ENDDO
     3045                ENDDO
     3046                sums_l(k,116,tn) = s1
    29833047             ENDDO
    29843048             !$acc end parallel loop
Note: See TracChangeset for help on using the changeset viewer.