Ignore:
Timestamp:
Oct 26, 2015 4:17:44 PM (8 years ago)
Author:
maronga
Message:

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

File:
1 edited

Legend:

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

    r1683 r1691  
    1414! PALM. If not, see <http://www.gnu.org/licenses/>.
    1515!
    16 ! Copyright 1997-2014 Leibniz Universitaet Hannover
     16! Copyright 1997-2015 Leibniz Universitaet Hannover
    1717!--------------------------------------------------------------------------------!
    1818!
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Revised calculation of Obukhov length. Added output of radiative heating >
     22! rates for RRTMG.
    2223!
    2324! Former revisions:
     
    161162
    162163    USE arrays_3d,                                                             &
    163         ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
    164                qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
    165                td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
    166                uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
     164        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, pt, q, qc, ql, qr, &
     165               qs, qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt,      &
     166               td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us,&
     167               usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
    167168       
    168169    USE cloud_parameters,                                                      &
     
    172173        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
    173174                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
    174                 large_scale_subsidence, max_pr_user, message_string, ocean,    &
    175                 passive_scalar, precipitation, simulated_time,                 &
     175                large_scale_subsidence, max_pr_user, message_string, neutral,  &
     176                ocean, passive_scalar, precipitation, simulated_time,          &
    176177                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
    177178                ws_scheme_mom, ws_scheme_sca
     
    199200    USE radiation_model_mod,                                                   &
    200201        ONLY:  dots_rad, radiation, radiation_scheme, rad_net,                 &
    201                rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
     202               rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr,                 &
     203               rad_sw_in, rad_sw_out, rad_sw_cs_hr, rad_sw_hr
    202204
    203205#if defined ( __rrtmg )
     
    207209 
    208210    USE statistics
     211
    209212
    210213    IMPLICIT NONE
     
    240243    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
    241244
    242     !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, usws, v, vsws, w )
     245    !$acc update host( km, kh, e, ol, pt, qs, qsws, shf, ts, u, usws, v, vsws, w )
    243246
    244247!
     
    744747!
    745748!--                   Formula does not work if ql(nzb) /= 0.0
    746                       sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
    747                                           qsws(j,i) * rmask(j,i,sr)
     749                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
     750                                          qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    748751                   ENDIF
    749752                ENDIF
    750753                IF ( passive_scalar )  THEN
    751                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
    752                                        qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     754                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
     755                                       qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    753756                ENDIF
    754757             ENDIF
     758
     759             IF ( .NOT. neutral )  THEN
     760                sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
     761                                    ol(j,i)  * rmask(j,i,sr) ! L
     762             ENDIF
     763
    755764
    756765             IF ( land_surface )  THEN
     
    774783#if defined ( __rrtmg )
    775784                IF ( radiation_scheme == 'rrtmg' )  THEN
    776                    sums_l(nzb,106,tn)  = sums_l(nzb,106,tn)  + rrtm_aldif(0,j,i)
    777                    sums_l(nzb,107,tn)  = sums_l(nzb,107,tn)  + rrtm_aldir(0,j,i)
    778                    sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  + rrtm_asdif(0,j,i)
    779                    sums_l(nzb,109,tn)  = sums_l(nzb,109,tn)  + rrtm_asdir(0,j,i)
     785                   sums_l(nzb,110,tn)  = sums_l(nzb,110,tn)  + rrtm_aldif(0,j,i)
     786                   sums_l(nzb,111,tn)  = sums_l(nzb,111,tn)  + rrtm_aldir(0,j,i)
     787                   sums_l(nzb,112,tn)  = sums_l(nzb,112,tn)  + rrtm_asdif(0,j,i)
     788                   sums_l(nzb,113,tn)  = sums_l(nzb,113,tn)  + rrtm_asdir(0,j,i)
    780789                ENDIF
    781790#endif
     
    9971006!--    First calculate the products, then the divergence.
    9981007!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
    999        IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
    1000 
     1008       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )   &
     1009       THEN
    10011010          sums_ll = 0.0_wp  ! local array
    10021011
     
    10371046!
    10381047!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
    1039        IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
    1040 
     1048       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )   &
     1049       THEN
    10411050          !$OMP DO
    10421051          DO  i = nxl, nxr
     
    11091118!--    Collect current large scale advection and subsidence tendencies for
    11101119!--    data output
    1111        IF ( large_scale_forcing  .AND.  ( simulated_time .GT. 0.0_wp ) )  THEN
     1120       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
    11121121!
    11131122!--       Interpolation in time of LSF_DATA
     
    11561165             DO  j =  nys, nyn
    11571166                DO  k = nzb_soil, nzt_soil
    1158                    sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i) * rmask(j,i,sr)
    1159                    sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i) * rmask(j,i,sr)
     1167                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
     1168                                      * rmask(j,i,sr)
     1169                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
     1170                                      * rmask(j,i,sr)
    11601171                ENDDO
    11611172             ENDDO
     
    11681179             DO  j =  nys, nyn
    11691180                DO  k = nzb_s_inner(j,i)+1, nzt+1
    1170                    sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i) * rmask(j,i,sr)
    1171                    sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i) * rmask(j,i,sr)
    1172                    sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i) * rmask(j,i,sr)
    1173                    sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i) * rmask(j,i,sr)
     1181                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
     1182                                       * rmask(j,i,sr)
     1183                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
     1184                                       * rmask(j,i,sr)
     1185                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
     1186                                       * rmask(j,i,sr)
     1187                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
     1188                                       * rmask(j,i,sr)
     1189                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
     1190                                       * rmask(j,i,sr)
     1191                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
     1192                                       * rmask(j,i,sr)
     1193                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
     1194                                       * rmask(j,i,sr)
     1195                   sums_l(k,109,tn)  = sums_l(k,108,tn)  + rad_sw_hr(k,j,i)    &
     1196                                       * rmask(j,i,sr)
    11741197                ENDDO
    11751198             ENDDO
     
    13741397          hom(:,1,105,sr) = sums(:,105)            ! rad_sw_out
    13751398
     1399          IF ( radiation_scheme == 'rrtmg' )  THEN
     1400             hom(:,1,106,sr) = sums(:,106)            ! rad_lw_cs_hr
     1401             hom(:,1,107,sr) = sums(:,107)            ! rad_lw_hr
     1402             hom(:,1,108,sr) = sums(:,108)            ! rad_sw_cs_hr
     1403             hom(:,1,109,sr) = sums(:,109)            ! rad_sw_hr
     1404
    13761405#if defined ( __rrtmg )
    1377           IF ( radiation_scheme == 'rrtmg' )  THEN
    1378              hom(:,1,106,sr) = sums(:,106)            ! rrtm_aldif
    1379              hom(:,1,107,sr) = sums(:,107)            ! rrtm_aldir
    1380              hom(:,1,108,sr) = sums(:,108)            ! rrtm_asdif
    1381              hom(:,1,109,sr) = sums(:,109)            ! rrtm_asdir
     1406             hom(:,1,110,sr) = sums(:,110)            ! rrtm_aldif
     1407             hom(:,1,111,sr) = sums(:,111)            ! rrtm_aldir
     1408             hom(:,1,112,sr) = sums(:,112)            ! rrtm_asdif
     1409             hom(:,1,113,sr) = sums(:,113)            ! rrtm_asdir
     1410#endif
    13821411          ENDIF
    1383 #endif
    1384 
    1385        ENDIF
     1412
     1413
     1414       ENDIF
     1415
     1416       hom(:,1,114,sr) = sums(:,114)            !: L
    13861417
    13871418       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
     
    15201551       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
    15211552       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
    1522 
     1553 
     1554!
     1555!--    Calculate obukhov length
    15231556       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
    1524           ts_value(22,sr) = ts_value(4,sr)**2 / &
    1525                             ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
     1557!           IF ( most_method == 'circular' )  THEN
     1558!              ts_value(22,sr) = ts_value(4,sr)**2 /                             &
     1559!                           ( kappa * g * ts_value(5,sr) / ts_value(18,sr) )
     1560!           ELSE
     1561             ts_value(22,sr) = hom(nzb,1,114,sr) 
     1562!          ENDIF     
    15261563       ELSE
    15271564          ts_value(22,sr) = 10000.0_wp
     
    15531590#if defined ( __rrtmg )
    15541591          IF ( radiation_scheme == 'rrtmg' )  THEN
    1555              ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
    1556              ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
    1557              ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
    1558              ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
     1592             ts_value(dots_rad+5,sr) = hom(nzb,1,110,sr)          ! rrtm_aldif
     1593             ts_value(dots_rad+6,sr) = hom(nzb,1,111,sr)          ! rrtm_aldir
     1594             ts_value(dots_rad+7,sr) = hom(nzb,1,112,sr)          ! rrtm_asdif
     1595             ts_value(dots_rad+8,sr) = hom(nzb,1,113,sr)          ! rrtm_asdir
    15591596          ENDIF
    15601597#endif
     
    16421679#if defined ( __rrtmg )
    16431680    USE radiation_model_mod,                                                   &
    1644         ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
     1681        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rad_lw_cs_hr,   &
     1682               rad_lw_hr,  rad_sw_cs_hr, rad_sw_hr
    16451683#endif
    16461684
     
    30633101!--    Collect current large scale advection and subsidence tendencies for
    30643102!--    data output
    3065        IF ( large_scale_forcing  .AND.  ( simulated_time .GT. 0.0_wp ) )  THEN
     3103       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
    30663104!
    30673105!--       Interpolation in time of LSF_DATA
     
    31103148             DO  j =  nys, nyn
    31113149                DO  k = nzb_soil, nzt_soil
    3112                    sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i) * rmask(j,i,sr)
    3113                    sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i) * rmask(j,i,sr)
     3150                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
     3151                                      * rmask(j,i,sr)
     3152                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
     3153                                      * rmask(j,i,sr)
    31143154                ENDDO
    31153155             ENDDO
     
    31233163             DO  j =  nys, nyn
    31243164                DO  k = nzb_s_inner(j,i)+1, nzt+1
    3125                    sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i) * rmask(j,i,sr)
    3126                    sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i) * rmask(j,i,sr)
    3127                    sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i) * rmask(j,i,sr)
    3128                    sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i) * rmask(j,i,sr)
     3165                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
     3166                                       * rmask(j,i,sr)
     3167                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
     3168                                       * rmask(j,i,sr)
     3169                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
     3170                                       * rmask(j,i,sr)
     3171                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
     3172                                       * rmask(j,i,sr)
     3173#if defined ( __rrtmg )
     3174                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
     3175                                       * rmask(j,i,sr)
     3176                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
     3177                                       * rmask(j,i,sr)
     3178                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
     3179                                       * rmask(j,i,sr)
     3180                   sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
     3181                                       * rmask(j,i,sr)
     3182#endif
    31293183                ENDDO
    31303184             ENDDO
     
    31933247          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
    31943248          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
    3195           sums(k,89:105)           = sums(k,89:105)     / ngp_2dh(sr)
    3196           sums(k,106:pr_palm-2)    = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
     3249          sums(k,89:114)           = sums(k,89:114)     / ngp_2dh(sr)
     3250          sums(k,115:pr_palm-2)    = sums(k,115:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
     3251
    31973252       ENDDO
    31983253
     
    34463501
    34473502       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
    3448           ts_value(22,sr) = ts_value(4,sr)**2_wp / &
    3449                             ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
     3503          ts_value(22,sr) = hom(nzb,1,114,sr)
    34503504       ELSE
    34513505          ts_value(22,sr) = 10000.0_wp
Note: See TracChangeset for help on using the changeset viewer.