Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 2032)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 2037)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Anelastic approximation implemented
 ! 
 ! Former revisions:
@@ -214,9 +214,11 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
-               qr, qs, qsws, qswst, rho_ocean, s, sa, ss, ssws, sswst, saswsb,       &
-               saswst, shf, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q,        &
-               time_vert, ts, tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, &
-               vswst, w, w_subs, zw
+        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
+               momentumflux_output_conversion, nr, ol, p, prho, prr, pt, q,    &
+               qc, ql, qr, qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, &
+               sa, ss, ssws, sswst, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q, &
+               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
+               uswst, vsws, v, vg, vpt, vswst, w, w_subs,                      &
+               waterflux_output_conversion, zw
         
     USE cloud_parameters,                                                      &
@@ -326,5 +328,6 @@
        sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
 !--    WARNING: next line still has to be adjusted for OpenMP 
-       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
+       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
+                        heatflux_output_conversion  ! heat flux from advec_s_bc
        sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
        sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
@@ -358,6 +361,8 @@
 !
 !--          Swap the turbulent quantities evaluated in advec_ws.
-             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u* 
-             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v* 
+             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
+                              * momentumflux_output_conversion ! w*u*
+             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
+                              * momentumflux_output_conversion ! w*v*
              sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2 
              sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2 
@@ -373,7 +378,9 @@
 
           DO  i = 0, threads_per_task-1
-             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i) ! w*pt* 
+             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)      &
+                                           * heatflux_output_conversion  ! w*pt*
              IF ( ocean          ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
-             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)  ! w*q*
+             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
+                                           * waterflux_output_conversion  ! w*q*
              IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
           ENDDO
@@ -722,5 +729,7 @@
                                    ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
                                  + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
-                                                               ) * rmask(j,i,sr)
+                                                           ) * rmask(j,i,sr)   &
+                                         * rho_air_zw(k)                       &
+                                         * momentumflux_output_conversion(k)
 !
 !--             Momentum flux w"v"
@@ -730,5 +739,7 @@
                                    ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
                                  + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
-                                                               ) * rmask(j,i,sr)
+                                                           ) * rmask(j,i,sr)   &
+                                         * rho_air_zw(k)                       &
+                                         * momentumflux_output_conversion(k)
 !
 !--             Heat flux w"pt"
@@ -736,4 +747,6 @@
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( pt(k+1,j,i) - pt(k,j,i) )   &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                * ddzu(k+1) * rmask(j,i,sr)
 
@@ -754,8 +767,12 @@
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                * ddzu(k+1) * rmask(j,i,sr)
                    sums_l(k,48,tn) = sums_l(k,48,tn)                           &
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
                                                * ( q(k+1,j,i) - q(k,j,i) )     &
+                                               * rho_air_zw(k)                 &
+                                               * waterflux_output_conversion(k)&
                                                * ddzu(k+1) * rmask(j,i,sr)
 
@@ -765,4 +782,6 @@
                                                * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
                                                 - ( q(k,j,i) - ql(k,j,i) ) )   &
+                                               * rho_air_zw(k)                 &
+                                               * waterflux_output_conversion(k)&
                                                * ddzu(k+1) * rmask(j,i,sr) 
                    ENDIF
@@ -784,8 +803,11 @@
              IF ( use_surface_fluxes )  THEN
                 sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
+                                    momentumflux_output_conversion(nzb) * &
                                     usws(j,i) * rmask(j,i,sr)     ! w"u"
                 sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
+                                    momentumflux_output_conversion(nzb) * &
                                     vsws(j,i) * rmask(j,i,sr)     ! w"v"
                 sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
+                                    heatflux_output_conversion(nzb) * &
                                     shf(j,i)  * rmask(j,i,sr)     ! w"pt"
                 sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
@@ -799,19 +821,23 @@
                 IF ( humidity )  THEN
                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
+                                       waterflux_output_conversion(nzb) *      &
                                        qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
                    sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
                                        ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
                                        shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
-                                                  qsws(j,i) )
+                                                  qsws(j,i) )                  &
+                                       * heatflux_output_conversion(nzb)
                    IF ( cloud_droplets )  THEN
                       sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
                                          ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
                                            ql(nzb,j,i) ) * shf(j,i) +          &
-                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) )
+                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) ) &
+                                          * heatflux_output_conversion(nzb)
                    ENDIF
                    IF ( cloud_physics )  THEN
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
-                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  & 
+                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
+                                          waterflux_output_conversion(nzb) *   &
                                           qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                    ENDIF
@@ -860,8 +886,11 @@
              IF ( use_top_fluxes )  THEN
                 sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
+                                    momentumflux_output_conversion(nzt:nzt+1) * &
                                     uswst(j,i) * rmask(j,i,sr)    ! w"u"
                 sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
+                                    momentumflux_output_conversion(nzt:nzt+1) * &
                                     vswst(j,i) * rmask(j,i,sr)    ! w"v"
                 sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
+                                    heatflux_output_conversion(nzt:nzt+1) * &
                                     tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
                 sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
@@ -876,14 +905,17 @@
                 IF ( humidity )  THEN
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
+                                       waterflux_output_conversion(nzt) *      &
                                        qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                    sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
                                        ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
                                        tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
-                                                             qswst(j,i) )
+                                                             qswst(j,i) )      &
+                                       * heatflux_output_conversion(nzt)
                    IF ( cloud_droplets )  THEN
                       sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
                                           ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
                                             ql(nzt,j,i) ) * tswst(j,i) +       &
-                                            0.61_wp * pt(nzt,j,i) * qswst(j,i) )
+                                           0.61_wp * pt(nzt,j,i) * qswst(j,i) )&
+                                           * heatflux_output_conversion(nzt)
                    ENDIF
                    IF ( cloud_physics )  THEN
@@ -891,4 +923,5 @@
 !--                   Formula does not work if ql(nzb) /= 0.0
                       sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
+                                          waterflux_output_conversion(nzt) *   &
                                           qswst(j,i) * rmask(j,i,sr)
                    ENDIF
@@ -943,4 +976,5 @@
                                     vpt(k+1,j,i) - hom(k+1,1,44,sr) )
                       sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
+                                               heatflux_output_conversion(k) * &
                                                           rmask(j,i,sr)
                       sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
@@ -953,4 +987,5 @@
                               hom(k+1,1,42,sr) )
                          sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
+                                             waterflux_output_conversion(k) *  &
                                                              rmask(j,i,sr)
                          sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
@@ -971,11 +1006,13 @@
                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) )
                          sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
+                                              heatflux_output_conversion(k) *  &
                                                              rmask(j,i,sr)
                       ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
-                         sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                & 
-                                             hom(k,1,41,sr) ) *                &
-                                           sums_l(k,17,tn) +                   &
-                                           0.61_wp * hom(k,1,4,sr) *           &
-                                           sums_l(k,49,tn)
+                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              & 
+                                               hom(k,1,41,sr) ) *              &
+                                             sums_l(k,17,tn) +                 &
+                                             0.61_wp * hom(k,1,4,sr) *         &
+                                             sums_l(k,49,tn)                   &
+                                           ) * heatflux_output_conversion(k)
                       END IF
                    END IF
@@ -996,4 +1033,5 @@
                 sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
                                              ( ust**2 + vst**2 + w(k,j,i)**2 ) &
+                                           * momentumflux_output_conversion(k) &
                                              * rmask(j,i,sr)
              ENDDO
@@ -1017,4 +1055,5 @@
                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
                                                     ( w(k,j,i-1) + w(k,j,i) )  &
+                                          * momentumflux_output_conversion(k)  &
                                                     * ust * rmask(j,i,sr)
 !
@@ -1022,4 +1061,5 @@
                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
                                                     ( w(k,j-1,i) + w(k,j,i) )  &
+                                          * momentumflux_output_conversion(k)  &
                                                     * vst * rmask(j,i,sr)
                ENDDO
@@ -1038,4 +1078,5 @@
                            ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
                              pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
+                           * heatflux_output_conversion(k)                     &
                            * w(k,j,i) * rmask(j,i,sr)
                   IF ( humidity )  THEN
@@ -1043,4 +1084,5 @@
                                       q(k+1,j,i) - hom(k+1,1,41,sr) )
                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
+                                       waterflux_output_conversion(k) *        &
                                        rmask(j,i,sr)
                   ENDIF
@@ -1147,8 +1189,12 @@
                                                    ( kh(k,j,i) + kh(k,j,i-1) ) &
                                                  * ( pt(k,j,i-1) - pt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                  * ddx * rmask(j,i,sr)
                    sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
                                                    ( kh(k,j,i) + kh(k,j-1,i) ) &
                                                  * ( pt(k,j-1,i) - pt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                  * ddy * rmask(j,i,sr)
 !
@@ -1157,5 +1203,6 @@
                                                   ( u(k,j,i) - hom(k,1,1,sr) ) &
                                     * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
-                                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
+                                               * heatflux_output_conversion(k)
                    pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
                                     pt(k,j,i)   - hom(k,1,4,sr) )
@@ -1163,5 +1210,6 @@
                                                   ( v(k,j,i) - hom(k,1,2,sr) ) &
                                     * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
-                                                 pt(k,j,i)   - hom(k,1,4,sr) )
+                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
+                                               * heatflux_output_conversion(k)
                 ENDDO
              ENDDO
@@ -1489,4 +1537,7 @@
        ENDIF
 
+       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
+       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
+
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
@@ -1597,5 +1648,6 @@
        THEN
           hom(nzb+8,1,pr_palm,sr) = &
-             ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
+             ( g / hom(k_surface_level+1,1,4,sr) *                             &
+             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
              * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
        ELSE
@@ -1708,9 +1760,10 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
-               qsws, qswst, rho_ocean, s, sa, saswsb, saswst, shf, ss, ssws, sswst,  &
-               td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts,      &
-               tswst, u, ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w,      &
-               w_subs, zw
+        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
+               momentumflux_output_conversion, nr, p, prho, pt, q, qc, ql, qr, &
+               qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, sa, saswsb, &
+               saswst, shf, ss, ssws, sswst, td_lsa_lpt, td_lsa_q, td_sub_lpt, &
+               td_sub_q, time_vert, ts, tswst, u, ug, us, usws, uswst, vsws,   &
+               v, vg, vpt, vswst, w, w_subs, waterflux_output_conversion, zw
                  
         
@@ -1828,5 +1881,6 @@
        sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
 !--    WARNING: next line still has to be adjusted for OpenMP
-       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
+       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
+                        heatflux_output_conversion  ! heat flux from advec_s_bc
        sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
        sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
@@ -1860,6 +1914,8 @@
 !
 !--          Swap the turbulent quantities evaluated in advec_ws.
-             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
-             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
+             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
+                              * momentumflux_output_conversion ! w*u*
+             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
+                              * momentumflux_output_conversion ! w*v*
              sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
              sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
@@ -1881,7 +1937,9 @@
 
           DO  i = 0, threads_per_task-1
-             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
+             sums_l(:,17,i) = sums_wspts_ws_l(:,i)                             &
+                              * heatflux_output_conversion        ! w*pt* from advec_s_ws
              IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
-             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i) !w*q*
+             IF ( humidity       )  sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)      &
+                                            * waterflux_output_conversion !w*q*
              IF ( passive_scalar )  sums_l(:,116,i) = sums_wsss_ws_l(:,i) !w*s*
           ENDDO
@@ -2378,5 +2436,7 @@
                                  + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
                                                                )               &
-                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
+                               * rho_air_zw(k)                                 &
+                               * momentumflux_output_conversion(k)
 !
 !--             Momentum flux w"v"
@@ -2387,9 +2447,13 @@
                                  + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
                                                                )               &
-                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)        &
+                               * rho_air_zw(k)                                 &
+                               * momentumflux_output_conversion(k)
 !
 !--             Heat flux w"pt"
                 s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
                                  * ( pt(k+1,j,i) - pt(k,j,i) )                 &
+                                 * rho_air_zw(k)                               &
+                                 * heatflux_output_conversion(k)               &
                                  * ddzu(k+1) * rmask(j,i,sr)                   &
                                  * rflags_invers(j,i,k+1)
@@ -2435,8 +2499,12 @@
                    s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
                                     * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
+                                    * rho_air_zw(k)                            &
+                                    * heatflux_output_conversion(k)            &
                                     * ddzu(k+1) * rmask(j,i,sr)                &
                                     * rflags_invers(j,i,k+1)
                    s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
                                     * ( q(k+1,j,i) - q(k,j,i) )                &
+                                    * rho_air_zw(k)                            &
+                                    * waterflux_output_conversion(k)           &
                                     * ddzu(k+1) * rmask(j,i,sr)                &
                                     * rflags_invers(j,i,k+1)
@@ -2459,4 +2527,6 @@
                                        *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
                                           - ( q(k,j,i) - ql(k,j,i) ) )         &
+                                       * rho_air_zw(k)                         &
+                                       * waterflux_output_conversion(k)        &
                                        * ddzu(k+1) * rmask(j,i,sr)             & 
                                        * rflags_invers(j,i,k+1)
@@ -2506,7 +2576,10 @@
 !
 !--             Subgridscale fluxes in the Prandtl layer
-                s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
-                s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
-                s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
+                s1 = s1 + usws(j,i) * momentumflux_output_conversion(nzb)      &
+                                    * rmask(j,i,sr) ! w"u"
+                s2 = s2 + vsws(j,i) * momentumflux_output_conversion(nzb)      &
+                                    * rmask(j,i,sr) ! w"v"
+                s3 = s3 + shf(j,i)  * heatflux_output_conversion(nzb)          &
+                                    * rmask(j,i,sr) ! w"pt"
                 s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
                 s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
@@ -2545,7 +2618,9 @@
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                   s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)      &
+                                       * rmask(j,i,sr) ! w"q" (w"qv")
                    s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
-                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
+                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )           &
+                             * heatflux_output_conversion(nzb)
                 ENDDO
              ENDDO
@@ -2564,5 +2639,6 @@
                       s1 = s1 + ( ( 1.0_wp +                                   &
                                     0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
-                                  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
+                                 shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )&
+                                * heatflux_output_conversion(nzb)
                    ENDDO
                 ENDDO
@@ -2582,5 +2658,6 @@
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
-                      s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
+                      s1 = s1 + qsws(j,i) * waterflux_output_conversion(nzb)   &
+                                          * rmask(j,i,sr)   ! w"q" (w"qv")
                    ENDDO
                 ENDDO
@@ -2624,7 +2701,10 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
-                s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
-                s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
+                s1 = s1 + uswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
+                                     * rmask(j,i,sr)    ! w"u"
+                s2 = s2 + vswst(j,i) * momentumflux_output_conversion(nzt:nzt+1) &
+                                     * rmask(j,i,sr)    ! w"v"
+                s3 = s3 + tswst(j,i) * heatflux_output_conversion(nzt:nzt+1)   &
+                                     * rmask(j,i,sr)    ! w"pt"
                 s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
                 s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
@@ -2663,7 +2743,9 @@
              DO  i = nxl, nxr
                 DO  j =  nys, nyn
-                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                   s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)     &
+                                        * rmask(j,i,sr) ! w"q" (w"qv")
                    s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
-                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )
+                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )          &
+                             * heatflux_output_conversion(nzt)
                 ENDDO
              ENDDO
@@ -2683,5 +2765,6 @@
                                     0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
                                   tswst(j,i) +                                 &
-                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )
+                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )         &
+                                * heatflux_output_conversion(nzt)
                    ENDDO
                 ENDDO
@@ -2701,5 +2784,6 @@
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
-                      s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                      s1 = s1 + qswst(j,i) * waterflux_output_conversion(nzt)  &
+                                           * rmask(j,i,sr)  ! w"q" (w"qv")
                    ENDDO
                 ENDDO
@@ -2755,5 +2839,6 @@
 !--             Energy flux w*e* (has to be adjusted?)
                 s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
-                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
+                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)    &
+                                   * momentumflux_output_conversion(k)
              ENDDO
           ENDDO
@@ -2822,4 +2907,5 @@
                       s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
                                            vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
+                                         heatflux_output_conversion(k) *       &
                                          w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                    ENDDO
@@ -2839,4 +2925,5 @@
                          s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
                                               ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
+                                          * waterflux_output_conversion(k)                      &
                                           * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                       ENDDO
@@ -2928,4 +3015,5 @@
                          s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
                                               vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
+                                          * heatflux_output_conversion(k)       &
                                           * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
                       ENDDO
@@ -2939,6 +3027,8 @@
                 !$acc parallel loop present( hom, sums_l )
                 DO  k = nzb, nzt_diff
-                   sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
-                                                0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
+                   sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * &
+                                       sums_l(k,17,tn) + 0.61_wp *             &
+                                       hom(k,1,4,sr) * sums_l(k,49,tn)         &
+                                     ) * heatflux_output_conversion(k)
                 ENDDO
                 !$acc end parallel loop
@@ -2993,4 +3083,5 @@
                    s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
                                     * ust * rmask(j,i,sr)                      &
+                                    * momentumflux_output_conversion(k)        &
                                     * rflags_invers(j,i,k+1)
 !
@@ -2998,9 +3089,10 @@
                    s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
                                     * vst * rmask(j,i,sr)                      & 
+                                    * momentumflux_output_conversion(k)        &
                                     * rflags_invers(j,i,k+1)
                 ENDDO
              ENDDO
              sums_l(k,13,tn) = s1
-             sums_l(k,15,tn) = s1
+             sums_l(k,15,tn) = s2
           ENDDO
           !$acc end parallel loop
@@ -3021,4 +3113,5 @@
                    s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
                                         pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
+                                    * heatflux_output_conversion(k)            &
                                     * w(k,j,i) * rmask(j,i,sr)                 & 
                                     * rflags_invers(j,i,k+1)
@@ -3039,4 +3132,5 @@
                       s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
                                            q(k+1,j,i) - hom(k+1,1,41,sr) )     &
+                                       * waterflux_output_conversion(k)        &
                                        * w(k,j,i) * rmask(j,i,sr)              &
                                        * rflags_invers(j,i,k+1)
@@ -3167,8 +3261,12 @@
                                                    ( kh(k,j,i) + kh(k,j,i-1) ) &
                                                  * ( pt(k,j,i-1) - pt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                  * ddx * rmask(j,i,sr)
                    sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
                                                    ( kh(k,j,i) + kh(k,j-1,i) ) &
                                                  * ( pt(k,j-1,i) - pt(k,j,i) ) &
+                                               * rho_air_zw(k)                 &
+                                               * heatflux_output_conversion(k) &
                                                  * ddy * rmask(j,i,sr)
 !
@@ -3177,5 +3275,6 @@
                                      ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
                                      ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
-                                       pt(k,j,i)   - hom(k,1,4,sr) )
+                                       pt(k,j,i)   - hom(k,1,4,sr) )           &
+                                     * heatflux_output_conversion(k)
                    pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
                                     pt(k,j,i)   - hom(k,1,4,sr) )
@@ -3183,5 +3282,6 @@
                                      ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
                                      ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
-                                       pt(k,j,i)   - hom(k,1,4,sr) )
+                                       pt(k,j,i)   - hom(k,1,4,sr) )           &
+                                     * heatflux_output_conversion(k)
                 ENDDO
              ENDDO
@@ -3468,4 +3568,7 @@
        END IF
 
+       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
+       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
+
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
@@ -3576,5 +3679,6 @@
        IF ( hom(nzb,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )  THEN
           hom(nzb+8,1,pr_palm,sr) = &
-             ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
+             ( g / hom(k_surface_level+1,1,4,sr) *                             &
+             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
              * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
        ELSE
