Ignore:
Timestamp:
May 30, 2017 5:47:52 PM (4 years ago)
Author:
suehring
Message:

Adjustments according new topography and surface-modelling concept implemented

File:
1 edited

Legend:

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

    r2119 r2232  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adjustments to new topography and surface concept
    2323!
    2424! Former revisions:
     
    224224    USE arrays_3d,                                                             &
    225225        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
    226                momentumflux_output_conversion, nr, ol, p, prho, prr, pt, q,    &
    227                qc, ql, qr, qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, &
    228                sa, ss, ssws, sswst, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q, &
    229                td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
    230                uswst, vsws, v, vg, vpt, vswst, w, w_subs,                      &
    231                waterflux_output_conversion, zw
     226               momentumflux_output_conversion, nr, p, prho, prr, pt, q,        &
     227               qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s,                  &
     228               sa, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, u,   &
     229               ug, v, vg, vpt, w, w_subs, waterflux_output_conversion, zw
    232230       
    233231    USE cloud_parameters,                                                      &
     
    236234    USE control_parameters,                                                    &
    237235        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    238                 dt_3d, g, humidity, kappa, large_scale_forcing,                &
     236                dt_3d, g, humidity, kappa, land_surface, large_scale_forcing,  &
    239237                large_scale_subsidence, max_pr_user, message_string, neutral,  &
    240238                microphysics_seifert, ocean, passive_scalar, simulated_time,   &
     
    250248    USE indices,                                                               &
    251249        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
    252                 ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
    253                 nzb_s_inner, nzt, nzt_diff
     250                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
    254251       
    255252    USE kinds
    256253   
    257254    USE land_surface_model_mod,                                                &
    258         ONLY:   ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,              &
    259                 qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
    260                 shf_eb, t_soil
     255        ONLY:   m_soil_h, nzb_soil, nzt_soil, t_soil_h
    261256
    262257    USE netcdf_interface,                                                      &
     
    277272    USE statistics
    278273
     274       USE surface_mod,                                                        &
     275          ONLY :  surf_def_h, surf_lsm_h, surf_usm_h
     276
    279277
    280278    IMPLICIT NONE
     
    283281    INTEGER(iwp) ::  j                   !<
    284282    INTEGER(iwp) ::  k                   !<
     283    INTEGER(iwp) ::  ki                  !<
    285284    INTEGER(iwp) ::  k_surface_level     !<
     285    INTEGER(iwp) ::  m                   !< loop variable over all horizontal wall elements
     286    INTEGER(iwp) ::  l                   !< loop variable over surface facing -- up- or downward-facing
    286287    INTEGER(iwp) ::  nt                  !<
    287288    INTEGER(iwp) ::  omp_get_thread_num  !<
    288289    INTEGER(iwp) ::  sr                  !<
    289290    INTEGER(iwp) ::  tn                  !<
    290    
     291
    291292    LOGICAL ::  first  !<
    292293   
    293294    REAL(wp) ::  dptdz_threshold  !<
    294295    REAL(wp) ::  fac              !<
     296    REAL(wp) ::  flag             !<
    295297    REAL(wp) ::  height           !<
    296298    REAL(wp) ::  pts              !<
     
    400402!--    for other horizontal averages.
    401403       tn = 0
    402 
    403        !$OMP PARALLEL PRIVATE( i, j, k, tn )
    404 !$     tn = omp_get_thread_num()
    405 
     404       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
     405       !$ tn = omp_get_thread_num()
    406406       !$OMP DO
    407407       DO  i = nxl, nxr
    408408          DO  j =  nys, nyn
    409              DO  k = nzb_s_inner(j,i), nzt+1
    410                 sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
    411                 sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
    412                 sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
     409             DO  k = nzb, nzt+1
     410                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     411                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)  &
     412                                                              * flag
     413                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)  &
     414                                                              * flag
     415                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)  &
     416                                                              * flag
    413417             ENDDO
    414418          ENDDO
     
    421425          DO  i = nxl, nxr
    422426             DO  j =  nys, nyn
    423                 DO  k = nzb_s_inner(j,i), nzt+1
    424                    sums_l(k,23,tn)  = sums_l(k,23,tn) + &
    425                                       sa(k,j,i) * rmask(j,i,sr)
     427                DO  k = nzb, nzt+1
     428                   sums_l(k,23,tn)  = sums_l(k,23,tn) + sa(k,j,i)              &
     429                                    * rmask(j,i,sr)                            &
     430                                    * MERGE( 1.0_wp, 0.0_wp,                   &
     431                                             BTEST( wall_flags_0(k,j,i), 22 ) )
    426432                ENDDO
    427433             ENDDO
     
    437443          DO  i = nxl, nxr
    438444             DO  j =  nys, nyn
    439                 DO  k = nzb_s_inner(j,i), nzt+1
    440                    sums_l(k,44,tn)  = sums_l(k,44,tn) + &
    441                                       vpt(k,j,i) * rmask(j,i,sr)
    442                    sums_l(k,41,tn)  = sums_l(k,41,tn) + &
    443                                       q(k,j,i) * rmask(j,i,sr)
     445                DO  k = nzb, nzt+1
     446                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     447                   sums_l(k,44,tn)  = sums_l(k,44,tn) +                        &
     448                                      vpt(k,j,i) * rmask(j,i,sr) * flag
     449                   sums_l(k,41,tn)  = sums_l(k,41,tn) +                        &
     450                                      q(k,j,i) * rmask(j,i,sr)   * flag
    444451                ENDDO
    445452             ENDDO
     
    449456             DO  i = nxl, nxr
    450457                DO  j =  nys, nyn
    451                    DO  k = nzb_s_inner(j,i), nzt+1
    452                       sums_l(k,42,tn) = sums_l(k,42,tn) + &
    453                                       ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
    454                       sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
     458                   DO  k = nzb, nzt+1
     459                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     460                      sums_l(k,42,tn) = sums_l(k,42,tn) +                      &
     461                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr) &
     462                                                               * flag
     463                      sums_l(k,43,tn) = sums_l(k,43,tn) + (                    &
    455464                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
    456                                                           ) * rmask(j,i,sr)
     465                                                          ) * rmask(j,i,sr)    &
     466                                                            * flag
    457467                   ENDDO
    458468                ENDDO
     
    467477          DO  i = nxl, nxr
    468478             DO  j =  nys, nyn
    469                 DO  k = nzb_s_inner(j,i), nzt+1
    470                    sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i) * rmask(j,i,sr)
     479                DO  k = nzb, nzt+1
     480                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i)             &
     481                                    * rmask(j,i,sr)                            &
     482                                    * MERGE( 1.0_wp, 0.0_wp,                   &
     483                                             BTEST( wall_flags_0(k,j,i), 22 ) )
    471484                ENDDO
    472485             ENDDO
     
    603616       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper,             &
    604617       !$OMP                   sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2,  &
    605        !$OMP                   w2 )
    606 !$     tn = omp_get_thread_num()
    607 
     618       !$OMP                   w2, flag, m, ki, l )
     619       !$ tn = omp_get_thread_num()
    608620       !$OMP DO
    609621       DO  i = nxl, nxr
    610622          DO  j =  nys, nyn
    611623             sums_l_etot = 0.0_wp
    612              DO  k = nzb_s_inner(j,i), nzt+1
     624             DO  k = nzb, nzt+1
     625                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
    613626!
    614627!--             Prognostic and diagnostic variables
    615                 sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
    616                 sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
    617                 sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
    618                 sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
    619                 sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
     628                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)  &
     629                                                              * flag
     630                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)  &
     631                                                              * flag
     632                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)  &
     633                                                              * flag
     634                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)  &
     635                                                              * flag
     636                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)  * flag
    620637
    621638                sums_l(k,33,tn) = sums_l(k,33,tn) + &
    622                                   ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
     639                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)&
     640                                                                 * flag
    623641
    624642                IF ( humidity )  THEN
    625643                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
    626                                   ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
     644                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)&
     645                                                                 * flag
    627646                ENDIF
    628647                IF ( passive_scalar )  THEN
    629648                   sums_l(k,118,tn) = sums_l(k,118,tn) + &
    630                                   ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)
     649                                  ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)&
     650                                                                  * flag
    631651                ENDIF
    632652!
    633653!--             Higher moments
    634654!--             (Computation of the skewness of w further below)
    635                 sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr)
     655                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) &
     656                                                                * flag
    636657
    637658                sums_l_etot  = sums_l_etot + &
    638                                         0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + &
    639                                         w(k,j,i)**2 ) * rmask(j,i,sr)
     659                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 +  &
     660                                        w(k,j,i)**2 )            * rmask(j,i,sr)&
     661                                                                 * flag
    640662             ENDDO
    641663!
     
    647669!
    648670!--          2D-arrays (being collected in the last column of sums_l)
    649              sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
    650                                         us(j,i)   * rmask(j,i,sr)
    651              sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
    652                                         usws(j,i) * rmask(j,i,sr)
    653              sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
    654                                         vsws(j,i) * rmask(j,i,sr)
    655              sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
    656                                         ts(j,i)   * rmask(j,i,sr)
    657              IF ( humidity )  THEN
    658                 sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
    659                                             qs(j,i)   * rmask(j,i,sr)
     671             IF ( surf_def_h(0)%ns >= 1 )  THEN
     672                m = surf_def_h(0)%start_index(j,i)
     673                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
     674                                        surf_def_h(0)%us(m)   * rmask(j,i,sr)
     675                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
     676                                        surf_def_h(0)%usws(m) * rmask(j,i,sr)
     677                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
     678                                        surf_def_h(0)%vsws(m) * rmask(j,i,sr)
     679                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
     680                                        surf_def_h(0)%ts(m)   * rmask(j,i,sr)
     681                IF ( humidity )  THEN
     682                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
     683                                            surf_def_h(0)%qs(m)   * rmask(j,i,sr)
     684                ENDIF
     685                IF ( passive_scalar )  THEN
     686                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
     687                                            surf_def_h(0)%ss(m)   * rmask(j,i,sr)
     688                ENDIF
    660689             ENDIF
    661              IF ( passive_scalar )  THEN
    662                 sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
    663                                             ss(j,i)   * rmask(j,i,sr)
     690             IF ( surf_lsm_h%ns >= 1 )  THEN
     691                m = surf_lsm_h%start_index(j,i)
     692                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
     693                                        surf_lsm_h%us(m)   * rmask(j,i,sr)
     694                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
     695                                        surf_lsm_h%usws(m) * rmask(j,i,sr)
     696                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
     697                                        surf_lsm_h%vsws(m) * rmask(j,i,sr)
     698                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
     699                                        surf_lsm_h%ts(m)   * rmask(j,i,sr)
     700                IF ( humidity )  THEN
     701                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
     702                                            surf_lsm_h%qs(m)   * rmask(j,i,sr)
     703                ENDIF
     704                IF ( passive_scalar )  THEN
     705                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
     706                                            surf_lsm_h%ss(m)   * rmask(j,i,sr)
     707                ENDIF
     708             ENDIF
     709             IF ( surf_usm_h%ns >= 1 )  THEN
     710                m = surf_lsm_h%start_index(j,i)
     711                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
     712                                        surf_usm_h%us(m)   * rmask(j,i,sr)
     713                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
     714                                        surf_usm_h%usws(m) * rmask(j,i,sr)
     715                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
     716                                        surf_usm_h%vsws(m) * rmask(j,i,sr)
     717                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
     718                                        surf_usm_h%ts(m)   * rmask(j,i,sr)
     719                IF ( humidity )  THEN
     720                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
     721                                            surf_usm_h%qs(m)   * rmask(j,i,sr)
     722                ENDIF
     723                IF ( passive_scalar )  THEN
     724                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
     725                                            surf_usm_h%ss(m)   * rmask(j,i,sr)
     726                ENDIF
    664727             ENDIF
    665728          ENDDO
     
    674737          DO  i = nxl, nxr
    675738             DO  j =  nys, nyn
    676                 DO  k = nzb_s_inner(j,i), nzt+1
     739                DO  k = nzb, nzt+1
     740                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     741
    677742                   u2   = u(k,j,i)**2
    678743                   v2   = v(k,j,i)**2
     
    681746                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
    682747
    683                    sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
    684                    sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
    685                    sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
     748                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)    &
     749                                                            * flag
     750                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)    &
     751                                                            * flag
     752                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)    &
     753                                                            * flag
    686754!
    687755!--                Perturbation energy
    688756
    689757                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
    690                                   ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
     758                                  ( ust2 + vst2 + w2 )      * rmask(j,i,sr)    &
     759                                                            * flag
    691760                ENDDO
    692761             ENDDO
     
    702771       DO  i = nxl, nxr
    703772          DO  j =  nys, nyn
    704              DO  k = nzb_s_inner(j,i), nzt+1
     773             DO  k = nzb, nzt+1
     774                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
     775
    705776                w2   = w(k,j,i)**2
    706777                ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
     
    709780
    710781                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
    711                                  + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
     782                                 + 0.5_wp * ( ust2 + vst2 + w2 )               &
     783                                 * rmask(j,i,sr)                               &
     784                                 * flag
    712785             ENDDO
    713786          ENDDO
     
    728801!--                However, this implies no error since staggered velocity
    729802!--                components are zero at the walls and inside buildings.
    730 
    731              DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
     803!--          Flag 23 is used to mask surface fluxes as well as model-top fluxes,
     804!--          which are added further below.
     805             DO  k = nzb, nzt
     806                flag = MERGE( 1.0_wp, 0.0_wp,                                  &
     807                              BTEST( wall_flags_0(k,j,i), 23 ) ) *             &
     808                       MERGE( 1.0_wp, 0.0_wp,                                  &
     809                              BTEST( wall_flags_0(k,j,i), 9  ) )
    732810!
    733811!--             Momentum flux w"u"
     
    739817                                                           ) * rmask(j,i,sr)   &
    740818                                         * rho_air_zw(k)                       &
    741                                          * momentumflux_output_conversion(k)
     819                                         * momentumflux_output_conversion(k)   &
     820                                         * flag
    742821!
    743822!--             Momentum flux w"v"
     
    749828                                                           ) * rmask(j,i,sr)   &
    750829                                         * rho_air_zw(k)                       &
    751                                          * momentumflux_output_conversion(k)
     830                                         * momentumflux_output_conversion(k)   &
     831                                         * flag
    752832!
    753833!--             Heat flux w"pt"
     
    757837                                               * rho_air_zw(k)                 &
    758838                                               * heatflux_output_conversion(k) &
    759                                                * ddzu(k+1) * rmask(j,i,sr)
     839                                               * ddzu(k+1) * rmask(j,i,sr)     &
     840                                               * flag
    760841
    761842
     
    766847                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    767848                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
    768                                                * ddzu(k+1) * rmask(j,i,sr)
     849                                               * ddzu(k+1) * rmask(j,i,sr)     &
     850                                               * flag
    769851                ENDIF
    770852
     
    777859                                               * rho_air_zw(k)                 &
    778860                                               * heatflux_output_conversion(k) &
    779                                                * ddzu(k+1) * rmask(j,i,sr)
     861                                               * ddzu(k+1) * rmask(j,i,sr) * flag
    780862                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
    781863                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
     
    783865                                               * rho_air_zw(k)                 &
    784866                                               * waterflux_output_conversion(k)&
    785                                                * ddzu(k+1) * rmask(j,i,sr)
     867                                               * ddzu(k+1) * rmask(j,i,sr) * flag
    786868
    787869                   IF ( cloud_physics ) THEN
     
    792874                                               * rho_air_zw(k)                 &
    793875                                               * waterflux_output_conversion(k)&
    794                                                * ddzu(k+1) * rmask(j,i,sr)
     876                                               * ddzu(k+1) * rmask(j,i,sr) * flag
    795877                   ENDIF
    796878                ENDIF
     
    802884                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    803885                                                  * ( s(k+1,j,i) - s(k,j,i) )  &
    804                                                   * ddzu(k+1) * rmask(j,i,sr)
     886                                                  * ddzu(k+1) * rmask(j,i,sr)  &
     887                                                              * flag
    805888                ENDIF
    806889
     
    810893!--          Subgridscale fluxes in the Prandtl layer
    811894             IF ( use_surface_fluxes )  THEN
    812                 sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
     895                DO  l = 0, 1
     896                   ki = MERGE( -1, 0, l == 0 )
     897                   IF ( surf_def_h(l)%ns >= 1 )  THEN
     898                      DO  m = surf_def_h(l)%start_index(j,i), surf_def_h(l)%end_index(j,i)
     899                         k = surf_def_h(l)%k(m)
     900
     901                         sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) + &
     902                                    momentumflux_output_conversion(k+ki) * &
     903                                    surf_def_h(l)%usws(m) * rmask(j,i,sr)     ! w"u"
     904                         sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) + &
     905                                    momentumflux_output_conversion(k+ki) * &
     906                                    surf_def_h(l)%vsws(m) * rmask(j,i,sr)     ! w"v"
     907                         sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) + &
     908                                    heatflux_output_conversion(k+ki) * &
     909                                    surf_def_h(l)%shf(m)  * rmask(j,i,sr)     ! w"pt"
     910                         sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) + &
     911                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
     912                         sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + &
     913                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
     914                         IF ( ocean )  THEN
     915                            sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + &
     916                                       surf_def_h(l)%sasws(m) * rmask(j,i,sr)  ! w"sa"
     917                         ENDIF
     918                         IF ( humidity )  THEN
     919                            sums_l(k+ki,48,tn) = sums_l(k+ki,48,tn) +                     &
     920                                       waterflux_output_conversion(k+ki) *      &
     921                                       surf_def_h(l)%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
     922                            sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                   &
     923                                       ( 1.0_wp + 0.61_wp * q(k+ki,j,i) ) *     &
     924                                       surf_def_h(l)%shf(m) + 0.61_wp * pt(k+ki,j,i) *      &
     925                                                  surf_def_h(l)%qsws(m) )                  &
     926                                       * heatflux_output_conversion(k+ki)
     927                            IF ( cloud_droplets )  THEN
     928                               sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                &
     929                                         ( 1.0_wp + 0.61_wp * q(k+ki,j,i) -     &
     930                                           ql(k+ki,j,i) ) * surf_def_h(l)%shf(m) +          &
     931                                           0.61_wp * pt(k+ki,j,i) * surf_def_h(l)%qsws(m) ) &
     932                                          * heatflux_output_conversion(k+ki)
     933                            ENDIF
     934                            IF ( cloud_physics )  THEN
     935!
     936!--                            Formula does not work if ql(k+ki) /= 0.0
     937                               sums_l(k+ki,51,tn) = sums_l(k+ki,51,tn) +                  &
     938                                          waterflux_output_conversion(k+ki) *   &
     939                                          surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
     940                            ENDIF
     941                         ENDIF
     942                         IF ( passive_scalar )  THEN
     943                            sums_l(k+ki,119,tn) = sums_l(k+ki,119,tn) +                     &
     944                                        surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
     945                         ENDIF
     946
     947                      ENDDO
     948
     949                   ENDIF
     950                ENDDO
     951                IF ( surf_lsm_h%ns >= 1 )  THEN
     952                   m = surf_lsm_h%start_index(j,i)
     953                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
    813954                                    momentumflux_output_conversion(nzb) * &
    814                                     usws(j,i) * rmask(j,i,sr)     ! w"u"
    815                 sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
     955                                    surf_lsm_h%usws(m) * rmask(j,i,sr)     ! w"u"
     956                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
    816957                                    momentumflux_output_conversion(nzb) * &
    817                                     vsws(j,i) * rmask(j,i,sr)     ! w"v"
    818                 sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
     958                                    surf_lsm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
     959                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
    819960                                    heatflux_output_conversion(nzb) * &
    820                                     shf(j,i)  * rmask(j,i,sr)     ! w"pt"
    821                 sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
     961                                    surf_lsm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
     962                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
    822963                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
    823                 sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
     964                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
    824965                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
    825                 IF ( ocean )  THEN
    826                    sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
    827                                        saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
    828                 ENDIF
    829                 IF ( humidity )  THEN
    830                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
     966                   IF ( ocean )  THEN
     967                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
     968                                       surf_lsm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
     969                   ENDIF
     970                   IF ( humidity )  THEN
     971                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
    831972                                       waterflux_output_conversion(nzb) *      &
    832                                        qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
    833                    sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
     973                                       surf_lsm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
     974                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
    834975                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
    835                                        shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
    836                                                   qsws(j,i) )                  &
     976                                       surf_lsm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
     977                                                  surf_lsm_h%qsws(m) )                  &
    837978                                       * heatflux_output_conversion(nzb)
    838                    IF ( cloud_droplets )  THEN
    839                       sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
     979                      IF ( cloud_droplets )  THEN
     980                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
    840981                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
    841                                            ql(nzb,j,i) ) * shf(j,i) +          &
    842                                            0.61_wp * pt(nzb,j,i) * qsws(j,i) ) &
     982                                           ql(nzb,j,i) ) * surf_lsm_h%shf(m) +          &
     983                                           0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) ) &
    843984                                          * heatflux_output_conversion(nzb)
     985                      ENDIF
     986                      IF ( cloud_physics )  THEN
     987!
     988!--                      Formula does not work if ql(nzb) /= 0.0
     989                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
     990                                          waterflux_output_conversion(nzb) *   &
     991                                          surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
     992                      ENDIF
    844993                   ENDIF
    845                    IF ( cloud_physics )  THEN
    846 !
    847 !--                   Formula does not work if ql(nzb) /= 0.0
    848                       sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
     994                   IF ( passive_scalar )  THEN
     995                      sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
     996                                        surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
     997                   ENDIF
     998
     999
     1000                ENDIF
     1001                IF ( surf_usm_h%ns >= 1 )  THEN
     1002                   m = surf_usm_h%start_index(j,i)
     1003                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
     1004                                    momentumflux_output_conversion(nzb) * &
     1005                                    surf_usm_h%usws(m) * rmask(j,i,sr)     ! w"u"
     1006                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
     1007                                    momentumflux_output_conversion(nzb) * &
     1008                                    surf_usm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
     1009                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
     1010                                    heatflux_output_conversion(nzb) * &
     1011                                    surf_usm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
     1012                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
     1013                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
     1014                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
     1015                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
     1016                   IF ( ocean )  THEN
     1017                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
     1018                                       surf_usm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
     1019                   ENDIF
     1020                   IF ( humidity )  THEN
     1021                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
     1022                                       waterflux_output_conversion(nzb) *      &
     1023                                       surf_usm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
     1024                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
     1025                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
     1026                                       surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
     1027                                                  surf_usm_h%qsws(m) )                  &
     1028                                       * heatflux_output_conversion(nzb)
     1029                      IF ( cloud_droplets )  THEN
     1030                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
     1031                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
     1032                                           ql(nzb,j,i) ) * surf_usm_h%shf(m) +          &
     1033                                           0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) ) &
     1034                                          * heatflux_output_conversion(nzb)
     1035                      ENDIF
     1036                      IF ( cloud_physics )  THEN
     1037!
     1038!--                      Formula does not work if ql(nzb) /= 0.0
     1039                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
    8491040                                          waterflux_output_conversion(nzb) *   &
    850                                           qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
     1041                                          surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
     1042                      ENDIF
    8511043                   ENDIF
    852                 ENDIF
    853                 IF ( passive_scalar )  THEN
    854                    sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
    855                                         ssws(j,i) * rmask(j,i,sr) ! w"s"
    856                 ENDIF
     1044                   IF ( passive_scalar )  THEN
     1045                      sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
     1046                                        surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s"
     1047                   ENDIF
     1048
     1049
     1050                ENDIF
     1051
    8571052             ENDIF
    8581053
    8591054             IF ( .NOT. neutral )  THEN
    860                 sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
    861                                     ol(j,i)  * rmask(j,i,sr) ! L
    862              ENDIF
    863 
    864 
    865              IF ( land_surface )  THEN
    866                 sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + ghf_eb(j,i)
    867                 sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + shf_eb(j,i)
    868                 sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + qsws_eb(j,i)
    869                 sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + qsws_liq_eb(j,i)
    870                 sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + qsws_soil_eb(j,i)
    871                 sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + qsws_veg_eb(j,i)
    872                 sums_l(nzb,99,tn)  = sums_l(nzb,99,tn) + r_a(j,i)
    873                 sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ r_s(j,i)
     1055                IF ( surf_def_h(0)%ns >= 1 )  THEN
     1056                   m = surf_def_h(0)%start_index(j,i)
     1057                   sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
     1058                                        surf_def_h(0)%ol(m)  * rmask(j,i,sr) ! L
     1059                ENDIF
     1060                IF ( surf_lsm_h%ns >= 1 )  THEN
     1061                   m = surf_lsm_h%start_index(j,i)
     1062                   sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
     1063                                        surf_lsm_h%ol(m)  * rmask(j,i,sr) ! L
     1064                ENDIF
     1065                IF ( surf_usm_h%ns >= 1 )  THEN
     1066                   m = surf_usm_h%start_index(j,i)
     1067                   sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
     1068                                        surf_usm_h%ol(m)  * rmask(j,i,sr) ! L
     1069                ENDIF
    8741070             ENDIF
    8751071
     
    8931089!--          Subgridscale fluxes at the top surface
    8941090             IF ( use_top_fluxes )  THEN
     1091                m = surf_def_h(2)%start_index(j,i)
    8951092                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
    8961093                                    momentumflux_output_conversion(nzt:nzt+1) * &
    897                                     uswst(j,i) * rmask(j,i,sr)    ! w"u"
     1094                                    surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
    8981095                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
    8991096                                    momentumflux_output_conversion(nzt:nzt+1) * &
    900                                     vswst(j,i) * rmask(j,i,sr)    ! w"v"
     1097                                    surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
    9011098                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
    9021099                                    heatflux_output_conversion(nzt:nzt+1) * &
    903                                     tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
     1100                                    surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
    9041101                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
    9051102                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
     
    9091106                IF ( ocean )  THEN
    9101107                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
    911                                        saswst(j,i) * rmask(j,i,sr)  ! w"sa"
     1108                                       surf_def_h(2)%sasws(m) * rmask(j,i,sr)  ! w"sa"
    9121109                ENDIF
    9131110                IF ( humidity )  THEN
    9141111                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
    9151112                                       waterflux_output_conversion(nzt) *      &
    916                                        qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
     1113                                       surf_def_h(2)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
    9171114                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
    9181115                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
    919                                        tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
    920                                                              qswst(j,i) )      &
     1116                                       surf_def_h(2)%shf(m) +                  &
     1117                                       0.61_wp * pt(nzt,j,i) *    &
     1118                                       surf_def_h(2)%qsws(m) )      &
    9211119                                       * heatflux_output_conversion(nzt)
    9221120                   IF ( cloud_droplets )  THEN
    9231121                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
    9241122                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
    925                                             ql(nzt,j,i) ) * tswst(j,i) +       &
    926                                            0.61_wp * pt(nzt,j,i) * qswst(j,i) )&
     1123                                            ql(nzt,j,i) ) *                    &
     1124                                            surf_def_h(2)%shf(m) +             &
     1125                                           0.61_wp * pt(nzt,j,i) *             &
     1126                                           surf_def_h(2)%qsws(m) )&
    9271127                                           * heatflux_output_conversion(nzt)
    9281128                   ENDIF
     
    9321132                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
    9331133                                          waterflux_output_conversion(nzt) *   &
    934                                           qswst(j,i) * rmask(j,i,sr)
     1134                                          surf_def_h(2)%qsws(m) * rmask(j,i,sr)
    9351135                   ENDIF
    9361136                ENDIF
    9371137                IF ( passive_scalar )  THEN
    9381138                   sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
    939                                         sswst(j,i) * rmask(j,i,sr) ! w"s"
     1139                                        surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s"
    9401140                ENDIF
    9411141             ENDIF
     
    9461146!--          ----  speaking the following k-loop would have to be split up and
    9471147!--                rearranged according to the staggered grid.
    948              DO  k = nzb_s_inner(j,i), nzt
     1148             DO  k = nzb, nzt
     1149                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
    9491150                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
    9501151                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
     
    9561157!--             Higher moments
    9571158                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
    958                                                     rmask(j,i,sr)
     1159                                                    rmask(j,i,sr) * flag
    9591160                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
    960                                                     rmask(j,i,sr)
     1161                                                    rmask(j,i,sr) * flag
    9611162
    9621163!
     
    9681169                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
    9691170                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
    970                                         rmask(j,i,sr)
     1171                                        rmask(j,i,sr) * flag
    9711172                   ENDIF
    972                    sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *            &
    973                                                        rmask(j,i,sr)
     1173                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *      &
     1174                                                       rmask(j,i,sr) * flag
    9741175                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
    975                                                        rmask(j,i,sr)
     1176                                                       rmask(j,i,sr) * flag
    9761177                ENDIF
    9771178
     
    9851186                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
    9861187                                               heatflux_output_conversion(k) * &
    987                                                           rmask(j,i,sr)
    988                       sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
     1188                                                          rmask(j,i,sr) * flag
     1189                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) &
     1190                                                                    * flag
    9891191
    9901192                      IF ( .NOT. cloud_droplets )  THEN
     
    9961198                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
    9971199                                             waterflux_output_conversion(k) *  &
    998                                                              rmask(j,i,sr)
     1200                                                             rmask(j,i,sr)  *  &
     1201                                                             flag
    9991202                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
    1000                                                              rmask(j,i,sr)
     1203                                                             rmask(j,i,sr) *   &
     1204                                                             flag
    10011205                         sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
    1002                                                              rmask(j,i,sr)
     1206                                                             rmask(j,i,sr) *   &
     1207                                                             flag
    10031208                         IF ( microphysics_seifert )  THEN
    10041209                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
    1005                                                                 rmask(j,i,sr)
     1210                                                                rmask(j,i,sr) *&
     1211                                                                flag
    10061212                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
    1007                                                                 rmask(j,i,sr)
     1213                                                                rmask(j,i,sr) *&
     1214                                                                flag
    10081215                         ENDIF
    10091216                      ENDIF
     
    10151222                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
    10161223                                              heatflux_output_conversion(k) *  &
    1017                                                              rmask(j,i,sr)
     1224                                                             rmask(j,i,sr)  *  &
     1225                                                             flag
    10181226                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
    10191227                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              &
     
    10221230                                             0.61_wp * hom(k,1,4,sr) *         &
    10231231                                             sums_l(k,49,tn)                   &
    1024                                            ) * heatflux_output_conversion(k)
     1232                                           ) * heatflux_output_conversion(k) * &
     1233                                               flag
    10251234                      END IF
    10261235                   END IF
     
    10331242                                    s(k+1,j,i) - hom(k+1,1,117,sr) )
    10341243                   sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
    1035                                                        rmask(j,i,sr)
     1244                                                       rmask(j,i,sr) * flag
    10361245                ENDIF
    10371246
     
    10421251                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
    10431252                                           * momentumflux_output_conversion(k) &
    1044                                              * rmask(j,i,sr)
     1253                                           * rmask(j,i,sr) * flag
    10451254             ENDDO
    10461255          ENDDO
    10471256       ENDDO
     1257       !$OMP END PARALLEL
     1258!
     1259!--    Treat land-surface quantities according to new wall model structure.
     1260       IF ( land_surface )  THEN
     1261          tn = 0
     1262          !$OMP PARALLEL PRIVATE( i, j, m, tn )
     1263          !$ tn = omp_get_thread_num()
     1264          !$OMP DO
     1265          DO  m = 1, surf_lsm_h%ns
     1266             i = surf_lsm_h%i(m)
     1267             j = surf_lsm_h%j(m)
     1268       
     1269             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
     1270                  j >= nys  .AND.  j <= nyn )  THEN
     1271                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf_eb(m)
     1272                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%shf_eb(m)
     1273                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + surf_lsm_h%qsws_eb(m)
     1274                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + surf_lsm_h%qsws_liq_eb(m)
     1275                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + surf_lsm_h%qsws_soil_eb(m)
     1276                sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + surf_lsm_h%qsws_veg_eb(m)
     1277                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn) + surf_lsm_h%r_a(m)
     1278                sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ surf_lsm_h%r_s(m)
     1279             ENDIF
     1280          ENDDO
     1281          !$OMP END PARALLEL
     1282
     1283          tn = 0
     1284          !$OMP PARALLEL PRIVATE( i, j, k, m, tn )
     1285          !$ tn = omp_get_thread_num()
     1286          !$OMP DO
     1287          DO  m = 1, surf_lsm_h%ns
     1288
     1289             i = surf_lsm_h%i(m)           
     1290             j = surf_lsm_h%j(m)
     1291
     1292             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
     1293                  j >= nys  .AND.  j <= nyn )  THEN
     1294
     1295                DO  k = nzb_soil, nzt_soil
     1296                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil_h%var_2d(k,m)  &
     1297                                      * rmask(j,i,sr)
     1298                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil_h%var_2d(k,m)  &
     1299                                      * rmask(j,i,sr)
     1300                ENDDO
     1301             ENDIF
     1302          ENDDO
     1303          !$OMP END PARALLEL
     1304       ENDIF
    10481305!
    10491306!--    For speed optimization fluxes which have been computed in part directly
    10501307!--    inside the WS advection routines are treated seperatly
    10511308!--    Momentum fluxes first:
     1309
     1310       tn = 0
     1311       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
     1312       !$ tn = omp_get_thread_num()
    10521313       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
    1053          !$OMP DO
    1054          DO  i = nxl, nxr
    1055             DO  j = nys, nyn
    1056                DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
    1057                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                &
    1058                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
    1059                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                &
    1060                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
    1061 !
    1062 !--               Momentum flux w*u*
    1063                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
    1064                                                     ( w(k,j,i-1) + w(k,j,i) )  &
    1065                                           * momentumflux_output_conversion(k)  &
    1066                                                     * ust * rmask(j,i,sr)
    1067 !
    1068 !--               Momentum flux w*v*
    1069                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
    1070                                                     ( w(k,j-1,i) + w(k,j,i) )  &
    1071                                           * momentumflux_output_conversion(k)  &
    1072                                                     * vst * rmask(j,i,sr)
    1073                ENDDO
    1074             ENDDO
    1075          ENDDO
     1314          !$OMP DO
     1315          DO  i = nxl, nxr
     1316             DO  j = nys, nyn
     1317                DO  k = nzb, nzt
     1318!
     1319!--                Flag 23 is used to mask surface fluxes as well as model-top
     1320!--                fluxes, which are added further below.
     1321                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
     1322                                 BTEST( wall_flags_0(k,j,i), 23 ) ) *          &
     1323                          MERGE( 1.0_wp, 0.0_wp,                               &
     1324                                 BTEST( wall_flags_0(k,j,i), 9  ) )
     1325
     1326                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
     1327                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
     1328                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
     1329                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
     1330!
     1331!--                Momentum flux w*u*
     1332                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                &
     1333                                                     ( w(k,j,i-1) + w(k,j,i) ) &
     1334                                           * momentumflux_output_conversion(k) &
     1335                                                     * ust * rmask(j,i,sr)     &
     1336                                                           * flag
     1337!
     1338!--                Momentum flux w*v*
     1339                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                &
     1340                                                     ( w(k,j-1,i) + w(k,j,i) ) &
     1341                                           * momentumflux_output_conversion(k) &
     1342                                                     * vst * rmask(j,i,sr)     &
     1343                                                           * flag
     1344                ENDDO
     1345             ENDDO
     1346          ENDDO
    10761347
    10771348       ENDIF
    10781349       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
    1079          !$OMP DO
    1080          DO  i = nxl, nxr
    1081             DO  j = nys, nyn
    1082                DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
    1083 !
    1084 !--               Vertical heat flux
    1085                   sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                 &
     1350          !$OMP DO
     1351          DO  i = nxl, nxr
     1352             DO  j = nys, nyn
     1353                DO  k = nzb, nzt
     1354                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
     1355                                 BTEST( wall_flags_0(k,j,i), 23 ) ) *          &
     1356                          MERGE( 1.0_wp, 0.0_wp,                               &
     1357                                 BTEST( wall_flags_0(k,j,i), 9  ) )
     1358!
     1359!--                Vertical heat flux
     1360                   sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                &
    10861361                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
    10871362                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
    10881363                           * heatflux_output_conversion(k)                     &
    1089                            * w(k,j,i) * rmask(j,i,sr)
    1090                   IF ( humidity )  THEN
    1091                      pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +            &
     1364                           * w(k,j,i) * rmask(j,i,sr) * flag
     1365                   IF ( humidity )  THEN
     1366                      pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +           &
    10921367                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
    1093                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
     1368                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *     &
    10941369                                       waterflux_output_conversion(k) *        &
    1095                                        rmask(j,i,sr)
    1096                   ENDIF
    1097                   IF ( passive_scalar )  THEN
    1098                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +            &
     1370                                       rmask(j,i,sr) * flag
     1371                   ENDIF
     1372                   IF ( passive_scalar )  THEN
     1373                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +           &
    10991374                                      s(k+1,j,i) - hom(k+1,1,117,sr) )
    1100                      sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *     &
    1101                                         rmask(j,i,sr)
    1102                   ENDIF
    1103                ENDDO
    1104             ENDDO
    1105          ENDDO
     1375                      sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *    &
     1376                                        rmask(j,i,sr) * flag
     1377                   ENDIF
     1378                ENDDO
     1379             ENDDO
     1380          ENDDO
    11061381
    11071382       ENDIF
     
    11261401          DO  i = nxl, nxr
    11271402             DO  j = nys, nyn
    1128                 DO  k = nzb_s_inner(j,i)+1, nzt
     1403                DO  k = nzb+1, nzt
     1404                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    11291405
    11301406                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
     
    11331409                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
    11341410                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
    1135                 + w(k,j,i)**2                                        )
     1411                + w(k,j,i)**2                                        ) * flag
    11361412
    11371413                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
    1138                                                * ( p(k,j,i) + p(k+1,j,i) )
     1414                                               * ( p(k,j,i) + p(k+1,j,i) ) * flag
    11391415
    11401416                ENDDO
     
    11641440          DO  i = nxl, nxr
    11651441             DO  j = nys, nyn
    1166                 DO  k = nzb_s_inner(j,i)+1, nzt
     1442                DO  k = nzb+1, nzt
     1443
     1444                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    11671445
    11681446                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
    11691447                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    11701448                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
    1171                                                                 ) * ddzw(k)
     1449                                                                ) * ddzw(k)    &
     1450                                                                  * flag
    11721451
    11731452                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
    11741453                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    1175                                                                 )
     1454                                                                )  * flag
    11761455
    11771456                ENDDO
     
    11911470          DO  i = nxl, nxr
    11921471             DO  j = nys, nyn
    1193                 DO  k = nzb_s_inner(j,i)+1, nzt
     1472                DO  k = nzb+1, nzt
     1473                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    11941474!
    11951475!--                Subgrid horizontal heat fluxes u"pt", v"pt"
     
    11991479                                               * rho_air_zw(k)                 &
    12001480                                               * heatflux_output_conversion(k) &
    1201                                                  * ddx * rmask(j,i,sr)
     1481                                                 * ddx * rmask(j,i,sr) * flag
    12021482                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
    12031483                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
     
    12051485                                               * rho_air_zw(k)                 &
    12061486                                               * heatflux_output_conversion(k) &
    1207                                                  * ddy * rmask(j,i,sr)
     1487                                                 * ddy * rmask(j,i,sr) * flag
    12081488!
    12091489!--                Resolved horizontal heat fluxes u*pt*, v*pt*
     
    12121492                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
    12131493                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
    1214                                                * heatflux_output_conversion(k)
     1494                                               * heatflux_output_conversion(k) &
     1495                                               * flag
    12151496                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
    12161497                                    pt(k,j,i)   - hom(k,1,4,sr) )
     
    12191500                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
    12201501                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
    1221                                                * heatflux_output_conversion(k)
     1502                                               * heatflux_output_conversion(k) &
     1503                                               * flag
    12221504                ENDDO
    12231505             ENDDO
     
    12791561       ENDIF
    12801562
     1563       tn = 0
    12811564       !$OMP PARALLEL PRIVATE( i, j, k, tn )
    1282 !$     tn = omp_get_thread_num()
    1283        IF ( land_surface )  THEN
    1284           !$OMP DO
    1285           DO  i = nxl, nxr
    1286              DO  j =  nys, nyn
    1287                 DO  k = nzb_soil, nzt_soil
    1288                    sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
    1289                                       * rmask(j,i,sr)
    1290                    sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
    1291                                       * rmask(j,i,sr)
    1292                 ENDDO
    1293              ENDDO
    1294           ENDDO
    1295        ENDIF
    1296        
     1565       !$ tn = omp_get_thread_num()       
    12971566       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
    12981567          !$OMP DO
    12991568          DO  i = nxl, nxr
    13001569             DO  j =  nys, nyn
    1301                 DO  k = nzb_s_inner(j,i)+1, nzt+1
     1570                DO  k = nzb+1, nzt+1
     1571                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     1572
    13021573                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
    1303                                        * rmask(j,i,sr)
     1574                                       * rmask(j,i,sr) * flag
    13041575                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
    1305                                        * rmask(j,i,sr)
     1576                                       * rmask(j,i,sr) * flag
    13061577                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
    1307                                        * rmask(j,i,sr)
     1578                                       * rmask(j,i,sr) * flag
    13081579                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
    1309                                        * rmask(j,i,sr)
     1580                                       * rmask(j,i,sr) * flag
    13101581                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
    1311                                        * rmask(j,i,sr)
     1582                                       * rmask(j,i,sr) * flag
    13121583                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
    1313                                        * rmask(j,i,sr)
     1584                                       * rmask(j,i,sr) * flag
    13141585                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
    1315                                        * rmask(j,i,sr)
     1586                                       * rmask(j,i,sr) * flag
    13161587                   sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
    1317                                        * rmask(j,i,sr)
     1588                                       * rmask(j,i,sr) * flag
    13181589                ENDDO
    13191590             ENDDO
Note: See TracChangeset for help on using the changeset viewer.