Ignore:
Timestamp:
Oct 26, 2016 11:15:40 AM (5 years ago)
Author:
knoop
Message:

Anelastic approximation implemented

File:
1 edited

Legend:

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

    r2032 r2037  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Anelastic approximation implemented
    2323!
    2424! Former revisions:
     
    214214
    215215    USE arrays_3d,                                                             &
    216         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
    217                qr, qs, qsws, qswst, rho_ocean, s, sa, ss, ssws, sswst, saswsb,       &
    218                saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q,        &
    219                time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, &
    220                vswst, w, w_subs, zw
     216        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
     217               momentumflux_output_conversion, nr, ol, p, prho, prr, pt, q,    &
     218               qc, ql, qr, qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, &
     219               sa, ss, ssws, sswst, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q, &
     220               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
     221               uswst, vsws, v, vg, vpt, vswst, w, w_subs,                      &
     222               waterflux_output_conversion, zw
    221223       
    222224    USE cloud_parameters,                                                      &
     
    326328       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
    327329!--    WARNING: next line still has to be adjusted for OpenMP
    328        sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
     330       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
     331                        heatflux_output_conversion  ! heat flux from advec_s_bc
    329332       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
    330333       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
     
    358361!
    359362!--          Swap the turbulent quantities evaluated in advec_ws.
    360              sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
    361              sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
     363             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
     364                              * momentumflux_output_conversion ! w*u*
     365             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
     366                              * momentumflux_output_conversion ! w*v*
    362367             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
    363368             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
     
    373378
    374379          DO  i = 0, threads_per_task-1
    375              sums_l(:,17,i)                        = sums_wspts_ws_l(:,i) ! w*pt*
     380             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)      &
     381                                           * heatflux_output_conversion  ! w*pt*
    376382             IF ( ocean          ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
    377              IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)  ! w*q*
     383             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
     384                                           * waterflux_output_conversion  ! w*q*
    378385             IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
    379386          ENDDO
     
    722729                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    723730                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
    724                                                                ) * rmask(j,i,sr)
     731                                                           ) * rmask(j,i,sr)   &
     732                                         * rho_air_zw(k)                       &
     733                                         * momentumflux_output_conversion(k)
    725734!
    726735!--             Momentum flux w"v"
     
    730739                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    731740                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    732                                                                ) * rmask(j,i,sr)
     741                                                           ) * rmask(j,i,sr)   &
     742                                         * rho_air_zw(k)                       &
     743                                         * momentumflux_output_conversion(k)
    733744!
    734745!--             Heat flux w"pt"
     
    736747                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    737748                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
     749                                               * rho_air_zw(k)                 &
     750                                               * heatflux_output_conversion(k) &
    738751                                               * ddzu(k+1) * rmask(j,i,sr)
    739752
     
    754767                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    755768                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
     769                                               * rho_air_zw(k)                 &
     770                                               * heatflux_output_conversion(k) &
    756771                                               * ddzu(k+1) * rmask(j,i,sr)
    757772                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
    758773                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    759774                                               * ( q(k+1,j,i) - q(k,j,i) )     &
     775                                               * rho_air_zw(k)                 &
     776                                               * waterflux_output_conversion(k)&
    760777                                               * ddzu(k+1) * rmask(j,i,sr)
    761778
     
    765782                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
    766783                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
     784                                               * rho_air_zw(k)                 &
     785                                               * waterflux_output_conversion(k)&
    767786                                               * ddzu(k+1) * rmask(j,i,sr)
    768787                   ENDIF
     
    784803             IF ( use_surface_fluxes )  THEN
    785804                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
     805                                    momentumflux_output_conversion(nzb) * &
    786806                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
    787807                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
     808                                    momentumflux_output_conversion(nzb) * &
    788809                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
    789810                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
     811                                    heatflux_output_conversion(nzb) * &
    790812                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
    791813                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
     
    799821                IF ( humidity )  THEN
    800822                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
     823                                       waterflux_output_conversion(nzb) *      &
    801824                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
    802825                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
    803826                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
    804827                                       shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
    805                                                   qsws(j,i) )
     828                                                  qsws(j,i) )                  &
     829                                       * heatflux_output_conversion(nzb)
    806830                   IF ( cloud_droplets )  THEN
    807831                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
    808832                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
    809833                                           ql(nzb,j,i) ) * shf(j,i) +          &
    810                                            0.61_wp * pt(nzb,j,i) * qsws(j,i) )
     834                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) ) &
     835                                          * heatflux_output_conversion(nzb)
    811836                   ENDIF
    812837                   IF ( cloud_physics )  THEN
    813838!
    814839!--                   Formula does not work if ql(nzb) /= 0.0
    815                       sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
     840                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
     841                                          waterflux_output_conversion(nzb) *   &
    816842                                          qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    817843                   ENDIF
     
    860886             IF ( use_top_fluxes )  THEN
    861887                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
     888                                    momentumflux_output_conversion(nzt:nzt+1) * &
    862889                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
    863890                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
     891                                    momentumflux_output_conversion(nzt:nzt+1) * &
    864892                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
    865893                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
     894                                    heatflux_output_conversion(nzt:nzt+1) * &
    866895                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
    867896                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
     
    876905                IF ( humidity )  THEN
    877906                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
     907                                       waterflux_output_conversion(nzt) *      &
    878908                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    879909                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
    880910                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
    881911                                       tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
    882                                                              qswst(j,i) )
     912                                                             qswst(j,i) )      &
     913                                       * heatflux_output_conversion(nzt)
    883914                   IF ( cloud_droplets )  THEN
    884915                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
    885916                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
    886917                                            ql(nzt,j,i) ) * tswst(j,i) +       &
    887                                             0.61_wp * pt(nzt,j,i) * qswst(j,i) )
     918                                           0.61_wp * pt(nzt,j,i) * qswst(j,i) )&
     919                                           * heatflux_output_conversion(nzt)
    888920                   ENDIF
    889921                   IF ( cloud_physics )  THEN
     
    891923!--                   Formula does not work if ql(nzb) /= 0.0
    892924                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
     925                                          waterflux_output_conversion(nzt) *   &
    893926                                          qswst(j,i) * rmask(j,i,sr)
    894927                   ENDIF
     
    943976                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    944977                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
     978                                               heatflux_output_conversion(k) * &
    945979                                                          rmask(j,i,sr)
    946980                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
     
    953987                              hom(k+1,1,42,sr) )
    954988                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
     989                                             waterflux_output_conversion(k) *  &
    955990                                                             rmask(j,i,sr)
    956991                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
     
    9711006                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    9721007                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
     1008                                              heatflux_output_conversion(k) *  &
    9731009                                                             rmask(j,i,sr)
    9741010                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
    975                          sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                &
    976                                              hom(k,1,41,sr) ) *                &
    977                                            sums_l(k,17,tn) +                   &
    978                                            0.61_wp * hom(k,1,4,sr) *           &
    979                                            sums_l(k,49,tn)
     1011                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              &
     1012                                               hom(k,1,41,sr) ) *              &
     1013                                             sums_l(k,17,tn) +                 &
     1014                                             0.61_wp * hom(k,1,4,sr) *         &
     1015                                             sums_l(k,49,tn)                   &
     1016                                           ) * heatflux_output_conversion(k)
    9801017                      END IF
    9811018                   END IF
     
    9961033                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
    9971034                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
     1035                                           * momentumflux_output_conversion(k) &
    9981036                                             * rmask(j,i,sr)
    9991037             ENDDO
     
    10171055                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
    10181056                                                    ( w(k,j,i-1) + w(k,j,i) )  &
     1057                                          * momentumflux_output_conversion(k)  &
    10191058                                                    * ust * rmask(j,i,sr)
    10201059!
     
    10221061                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
    10231062                                                    ( w(k,j-1,i) + w(k,j,i) )  &
     1063                                          * momentumflux_output_conversion(k)  &
    10241064                                                    * vst * rmask(j,i,sr)
    10251065               ENDDO
     
    10381078                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
    10391079                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
     1080                           * heatflux_output_conversion(k)                     &
    10401081                           * w(k,j,i) * rmask(j,i,sr)
    10411082                  IF ( humidity )  THEN
     
    10431084                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
    10441085                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
     1086                                       waterflux_output_conversion(k) *        &
    10451087                                       rmask(j,i,sr)
    10461088                  ENDIF
     
    11471189                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
    11481190                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
     1191                                               * rho_air_zw(k)                 &
     1192                                               * heatflux_output_conversion(k) &
    11491193                                                 * ddx * rmask(j,i,sr)
    11501194                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
    11511195                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
    11521196                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
     1197                                               * rho_air_zw(k)                 &
     1198                                               * heatflux_output_conversion(k) &
    11531199                                                 * ddy * rmask(j,i,sr)
    11541200!
     
    11571203                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
    11581204                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
    1159                                                  pt(k,j,i)   - hom(k,1,4,sr) )
     1205                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
     1206                                               * heatflux_output_conversion(k)
    11601207                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
    11611208                                    pt(k,j,i)   - hom(k,1,4,sr) )
     
    11631210                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
    11641211                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
    1165                                                  pt(k,j,i)   - hom(k,1,4,sr) )
     1212                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
     1213                                               * heatflux_output_conversion(k)
    11661214                ENDDO
    11671215             ENDDO
     
    14891537       ENDIF
    14901538
     1539       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
     1540       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
     1541
    14911542       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
    14921543                                       ! u*, w'u', w'v', t* (in last profile)
     
    15971648       THEN
    15981649          hom(nzb+8,1,pr_palm,sr) = &
    1599              ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
     1650             ( g / hom(k_surface_level+1,1,4,sr) *                             &
     1651             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
    16001652             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
    16011653       ELSE
     
    17081760
    17091761    USE arrays_3d,                                                             &
    1710         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
    1711                qsws, qswst, rho_ocean, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
    1712                td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts,      &
    1713                tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w,      &
    1714                w_subs, zw
     1762        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
     1763               momentumflux_output_conversion, nr, p, prho, pt, q, qc, ql, qr, &
     1764               qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, sa, saswsb, &
     1765               saswst, shf, ss, ssws, sswst, td_lsa_lpt, td_lsa_q, td_sub_lpt, &
     1766               td_sub_q, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,   &
     1767               v, vg, vpt, vswst, w, w_subs, waterflux_output_conversion, zw
    17151768                 
    17161769       
     
    18281881       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
    18291882!--    WARNING: next line still has to be adjusted for OpenMP
    1830        sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
     1883       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
     1884                        heatflux_output_conversion  ! heat flux from advec_s_bc
    18311885       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
    18321886       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
     
    18601914!
    18611915!--          Swap the turbulent quantities evaluated in advec_ws.
    1862              sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
    1863              sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
     1916             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
     1917                              * momentumflux_output_conversion ! w*u*
     1918             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
     1919                              * momentumflux_output_conversion ! w*v*
    18641920             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
    18651921             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
     
    18811937
    18821938          DO  i = 0, threads_per_task-1
    1883              sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
     1939             sums_l(:,17,i) = sums_wspts_ws_l(:,i)                             &
     1940                              * heatflux_output_conversion        ! w*pt* from advec_s_ws
    18841941             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
    1885              IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i) !w*q*
     1942             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)      &
     1943                                            * waterflux_output_conversion !w*q*
    18861944             IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
    18871945          ENDDO
     
    23782436                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
    23792437                                                               )               &
    2380                                * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2438                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
     2439                               * rho_air_zw(k)                                 &
     2440                               * momentumflux_output_conversion(k)
    23812441!
    23822442!--             Momentum flux w"v"
     
    23872447                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    23882448                                                               )               &
    2389                                * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2449                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
     2450                               * rho_air_zw(k)                                 &
     2451                               * momentumflux_output_conversion(k)
    23902452!
    23912453!--             Heat flux w"pt"
    23922454                s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
    23932455                                 * ( pt(k+1,j,i) - pt(k,j,i) )                 &
     2456                                 * rho_air_zw(k)                               &
     2457                                 * heatflux_output_conversion(k)               &
    23942458                                 * ddzu(k+1) * rmask(j,i,sr)                   &
    23952459                                 * rflags_invers(j,i,k+1)
     
    24352499                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    24362500                                    * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
     2501                                    * rho_air_zw(k)                            &
     2502                                    * heatflux_output_conversion(k)            &
    24372503                                    * ddzu(k+1) * rmask(j,i,sr)                &
    24382504                                    * rflags_invers(j,i,k+1)
    24392505                   s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
    24402506                                    * ( q(k+1,j,i) - q(k,j,i) )                &
     2507                                    * rho_air_zw(k)                            &
     2508                                    * waterflux_output_conversion(k)           &
    24412509                                    * ddzu(k+1) * rmask(j,i,sr)                &
    24422510                                    * rflags_invers(j,i,k+1)
     
    24592527                                       *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
    24602528                                          - ( q(k,j,i) - ql(k,j,i) ) )         &
     2529                                       * rho_air_zw(k)                         &
     2530                                       * waterflux_output_conversion(k)        &
    24612531                                       * ddzu(k+1) * rmask(j,i,sr)             &
    24622532                                       * rflags_invers(j,i,k+1)
     
    25062576!
    25072577!--             Subgridscale fluxes in the Prandtl layer
    2508                 s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
    2509                 s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
    2510                 s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
     2578                s1 = s1 + usws(j,i) * momentumflux_output_conversion(nzb)      &
     2579                                    * rmask(j,i,sr) ! w"u"
     2580                s2 = s2 + vsws(j,i) * momentumflux_output_conversion(nzb)      &
     2581                                    * rmask(j,i,sr) ! w"v"
     2582                s3 = s3 + shf(j,i)  * heatflux_output_conversion(nzb)          &
     2583                                    * rmask(j,i,sr) ! w"pt"
    25112584                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
    25122585                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
     
    25452618             DO  i = nxl, nxr
    25462619                DO  j =  nys, nyn
    2547                    s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     2620                   s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)      &
     2621                                       * rmask(j,i,sr) ! w"q" (w"qv")
    25482622                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
    2549                                + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
     2623                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )           &
     2624                             * heatflux_output_conversion(nzb)
    25502625                ENDDO
    25512626             ENDDO
     
    25642639                      s1 = s1 + ( ( 1.0_wp +                                   &
    25652640                                    0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
    2566                                   shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
     2641                                 shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )&
     2642                                * heatflux_output_conversion(nzb)
    25672643                   ENDDO
    25682644                ENDDO
     
    25822658!
    25832659!--                   Formula does not work if ql(nzb) /= 0.0
    2584                       s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
     2660                      s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)   &
     2661                                          * rmask(j,i,sr)   ! w"q" (w"qv")
    25852662                   ENDDO
    25862663                ENDDO
     
    26242701          DO  i = nxl, nxr
    26252702             DO  j =  nys, nyn
    2626                 s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
    2627                 s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
    2628                 s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
     2703                s1 = s1 + uswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
     2704                                     * rmask(j,i,sr)    ! w"u"
     2705                s2 = s2 + vswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
     2706                                     * rmask(j,i,sr)    ! w"v"
     2707                s3 = s3 + tswst(j,i) * heatflux_output_conversion(nzt:nzt+1)   &
     2708                                     * rmask(j,i,sr)    ! w"pt"
    26292709                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
    26302710                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
     
    26632743             DO  i = nxl, nxr
    26642744                DO  j =  nys, nyn
    2665                    s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
     2745                   s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)     &
     2746                                        * rmask(j,i,sr) ! w"q" (w"qv")
    26662747                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
    2667                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )
     2748                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )          &
     2749                             * heatflux_output_conversion(nzt)
    26682750                ENDDO
    26692751             ENDDO
     
    26832765                                    0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
    26842766                                  tswst(j,i) +                                 &
    2685                                   0.61_wp * pt(nzt,j,i) * qswst(j,i) )
     2767                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )         &
     2768                                * heatflux_output_conversion(nzt)
    26862769                   ENDDO
    26872770                ENDDO
     
    27012784!
    27022785!--                   Formula does not work if ql(nzb) /= 0.0
    2703                       s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     2786                      s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)  &
     2787                                           * rmask(j,i,sr)  ! w"q" (w"qv")
    27042788                   ENDDO
    27052789                ENDDO
     
    27552839!--             Energy flux w*e* (has to be adjusted?)
    27562840                s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
    2757                                    * rmask(j,i,sr) * rflags_invers(j,i,k+1)
     2841                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)    &
     2842                                   * momentumflux_output_conversion(k)
    27582843             ENDDO
    27592844          ENDDO
     
    28222907                      s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
    28232908                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
     2909                                         heatflux_output_conversion(k) *       &
    28242910                                         w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    28252911                   ENDDO
     
    28392925                         s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
    28402926                                              ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
     2927                                          * waterflux_output_conversion(k)                      &
    28412928                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    28422929                      ENDDO
     
    29283015                         s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
    29293016                                              vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
     3017                                          * heatflux_output_conversion(k)       &
    29303018                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
    29313019                      ENDDO
     
    29393027                !$acc parallel loop present( hom, sums_l )
    29403028                DO  k = nzb, nzt_diff
    2941                    sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
    2942                                                 0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
     3029                   sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * &
     3030                                       sums_l(k,17,tn) + 0.61_wp *             &
     3031                                       hom(k,1,4,sr) * sums_l(k,49,tn)         &
     3032                                     ) * heatflux_output_conversion(k)
    29433033                ENDDO
    29443034                !$acc end parallel loop
     
    29933083                   s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
    29943084                                    * ust * rmask(j,i,sr)                      &
     3085                                    * momentumflux_output_conversion(k)        &
    29953086                                    * rflags_invers(j,i,k+1)
    29963087!
     
    29983089                   s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
    29993090                                    * vst * rmask(j,i,sr)                      &
     3091                                    * momentumflux_output_conversion(k)        &
    30003092                                    * rflags_invers(j,i,k+1)
    30013093                ENDDO
    30023094             ENDDO
    30033095             sums_l(k,13,tn) = s1
    3004              sums_l(k,15,tn) = s1
     3096             sums_l(k,15,tn) = s2
    30053097          ENDDO
    30063098          !$acc end parallel loop
     
    30213113                   s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
    30223114                                        pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
     3115                                    * heatflux_output_conversion(k)            &
    30233116                                    * w(k,j,i) * rmask(j,i,sr)                 &
    30243117                                    * rflags_invers(j,i,k+1)
     
    30393132                      s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
    30403133                                           q(k+1,j,i) - hom(k+1,1,41,sr) )     &
     3134                                       * waterflux_output_conversion(k)        &
    30413135                                       * w(k,j,i) * rmask(j,i,sr)              &
    30423136                                       * rflags_invers(j,i,k+1)
     
    31673261                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
    31683262                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
     3263                                               * rho_air_zw(k)                 &
     3264                                               * heatflux_output_conversion(k) &
    31693265                                                 * ddx * rmask(j,i,sr)
    31703266                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
    31713267                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
    31723268                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
     3269                                               * rho_air_zw(k)                 &
     3270                                               * heatflux_output_conversion(k) &
    31733271                                                 * ddy * rmask(j,i,sr)
    31743272!
     
    31773275                                     ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
    31783276                                     ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
    3179                                        pt(k,j,i)   - hom(k,1,4,sr) )
     3277                                       pt(k,j,i)   - hom(k,1,4,sr) )           &
     3278                                     * heatflux_output_conversion(k)
    31803279                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
    31813280                                    pt(k,j,i)   - hom(k,1,4,sr) )
     
    31833282                                     ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
    31843283                                     ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
    3185                                        pt(k,j,i)   - hom(k,1,4,sr) )
     3284                                       pt(k,j,i)   - hom(k,1,4,sr) )           &
     3285                                     * heatflux_output_conversion(k)
    31863286                ENDDO
    31873287             ENDDO
     
    34683568       END IF
    34693569
     3570       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
     3571       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
     3572
    34703573       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
    34713574                                       ! u*, w'u', w'v', t* (in last profile)
     
    35763679       IF ( hom(nzb,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )  THEN
    35773680          hom(nzb+8,1,pr_palm,sr) = &
    3578              ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
     3681             ( g / hom(k_surface_level+1,1,4,sr) *                             &
     3682             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
    35793683             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
    35803684       ELSE
Note: See TracChangeset for help on using the changeset viewer.