Changeset 1709


Ignore:
Timestamp:
Nov 4, 2015 2:47:01 PM (8 years ago)
Author:
maronga
Message:

several bugfixes related to the new surface layer routine and land-surface-radiation interaction

Location:
palm/trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/DOC/palm2doxygen.config

    r1684 r1709  
    812812# exclude all test directories for example use the pattern */test/*
    813813
    814 EXCLUDE_PATTERNS       = */message.f90 */cpulog.f90
     814EXCLUDE_PATTERNS       = */message.f90 */cpulog.f90  */tridia_solver.f90
    815815
    816816# The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names
  • palm/trunk/SOURCE/flow_statistics.f90

    r1701 r1709  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Updated output of Obukhov length
    2222!
    2323! Former revisions:
     
    15541554       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
    15551555       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
    1556  
    1557 !
    1558 !--    Calculate obukhov length
    1559        IF ( ts_value(5,sr) /= 0.0_wp )  THEN
    1560 !           IF ( most_method == 'circular' )  THEN
    1561 !              ts_value(22,sr) = ts_value(4,sr)**2 /                             &
    1562 !                           ( kappa * g * ts_value(5,sr) / ts_value(18,sr) )
    1563 !           ELSE
    1564              ts_value(22,sr) = hom(nzb,1,114,sr) 
    1565 !          ENDIF     
     1556
     1557       IF ( .NOT. neutral )  THEN
     1558          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
    15661559       ELSE
    1567           ts_value(22,sr) = 10000.0_wp
     1560          ts_value(22,sr) = 1.0E10_wp
    15681561       ENDIF
    15691562
     
    35033496       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
    35043497
    3505        IF ( ts_value(5,sr) /= 0.0_wp )  THEN
    3506           ts_value(22,sr) = hom(nzb,1,114,sr)
     3498       IF ( .NOT. neutral )  THEN
     3499          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
    35073500       ELSE
    3508           ts_value(22,sr) = 10000.0_wp
     3501          ts_value(22,sr) = 1.0E10_wp
    35093502       ENDIF
    35103503
  • palm/trunk/SOURCE/init_1d_model.f90

    r1698 r1709  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Set initial time step to 10 s to avoid instability of the 1d model for small
     22! grid spacings
    2223!
    2324! Former revisions:
     
    8182!>
    8283!> @todo harmonize code with new surface_layer_fluxes module
     84!> @bug 1D model crashes when using small grid spacings in the order of 1 m
    8385!------------------------------------------------------------------------------!
    8486 SUBROUTINE init_1d_model
     
    262264!-- Determine the time step at the start of a 1D-simulation and
    263265!-- determine and printout quantities used for run control
    264     CALL timestep_1d
     266    dt_1d = 10.0_wp
    265267    CALL run_control_1d
    266268
     
    864866    USE pegrid
    865867   
    866     USE control_parameters,                                                              &
     868    USE control_parameters,                                                    &
    867869        ONLY:  message_string
    868870
  • palm/trunk/SOURCE/land_surface_model.f90

    r1698 r1709  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Renamed pt_1 and qv_1 to pt1 and qv1.
     22! Bugfix: set initial values for t_surface_p in case of restart runs
     23! Bugfix: zero resistance caused crash when using radiation_scheme = 'clear-sky'
     24! Bugfix: calculation of rad_net when using radiation_scheme = 'clear-sky'
     25! Added todo action
    2226!
    2327! Former revisions:
     
    8791!>
    8892!> @todo Consider partial absorption of the net shortwave radiation by the
    89 !>       surface layer.
     93!>       skin layer.
    9094!> @todo Allow for water surfaces
    9195!> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0,
     
    9397!> @todo Implement surface runoff model (required when performing long-term LES
    9498!>       with considerable precipitation.
     99!> @todo Fix crashes with radiation_scheme == 'constant'
    95100!>
    96 !> @note  No time step criterion is required as long as the soil layers do not
    97 !>        become too thin.
     101!> @note No time step criterion is required as long as the soil layers do not
     102!>       become too thin.
    98103!------------------------------------------------------------------------------!
    99104 MODULE land_surface_model_mod
     
    124129    USE radiation_model_mod,                                                   &
    125130        ONLY:  force_radiation_call, radiation_scheme, rad_net, rad_sw_in,     &
    126                rad_lw_out, sigma_sb
     131               rad_lw_out, rad_lw_out_change_0, sigma_sb
    127132       
    128133#if defined ( __rrtmg )
    129134    USE radiation_model_mod,                                                   &
    130         ONLY:  rrtm_lwuflx_dt, rrtm_idrv
     135        ONLY:  rrtm_idrv
    131136#endif
    132137
     
    275280              rad_net_l,        & !< local copy of rad_net (net radiation at surface)
    276281              r_a,              & !< aerodynamic resistance
    277               r_a_av,           & !< avergae of r_a
     282              r_a_av,           & !< average of r_a
    278283              r_canopy,         & !< canopy resistance
    279284              r_soil,           & !< soil resistance
     
    635640       INTEGER(iwp) ::  k !< running index
    636641
    637        REAL(wp) :: pt_1   !< potential temperature at first grid level
     642       REAL(wp) :: pt1   !< potential temperature at first grid level
    638643
    639644
     
    788793!--       Calculate surface temperature
    789794          t_surface   = pt_surface * exn
    790           t_surface_p = t_surface
    791795
    792796!
     
    798802
    799803                IF ( cloud_physics )  THEN
    800                    pt_1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
     804                   pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
    801805                ELSE
    802                    pt_1 = pt(k+1,j,i)
     806                   pt1 = pt(k+1,j,i)
    803807                ENDIF
    804808
    805809!
    806810!--             Assure that r_a cannot be zero at model start
    807                 IF ( pt_1 == pt(k,j,i) )  pt_1 = pt_1 + 1.0E-10_wp
     811                IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
    808812
    809813                us(j,i)  = 0.1_wp
    810                 ts(j,i)  = (pt_1 - pt(k,j,i)) / r_a(j,i)
     814                ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
    811815                shf(j,i) = - us(j,i) * ts(j,i)
    812816             ENDDO
     
    899903       m_soil_p    = m_soil
    900904       m_liq_eb_p  = m_liq_eb
     905       t_surface_p = t_surface
     906
     907
    901908
    902909!--    Store initial profiles of t_soil and m_soil (assuming they are
     
    926933       dots_unit(dots_soil+6:dots_soil+7) = "s/m"
    927934
    928        RETURN
    929935
    930936    END SUBROUTINE init_lsm
     
    964970                   lambda_surface, & !< Current value of lambda_surface
    965971                   m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
    966                    pt_1,        & !< potential temperature at first grid level
    967                    qv_1           !< specific humidity at first grid level
     972                   pt1,         & !< potential temperature at first grid level
     973                   qv1            !< specific humidity at first grid level
    968974
    969975!
     
    989995!--          equivalent to the ECMWF formulation using drag coefficients
    990996             IF ( cloud_physics )  THEN
    991                 pt_1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
    992                 qv_1 = q(k+1,j,i) - ql(k+1,j,i)
     997                pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
     998                qv1 = q(k+1,j,i) - ql(k+1,j,i)
    993999             ELSE
    994                 pt_1 = pt(k+1,j,i)
    995                 qv_1 = q(k+1,j,i)
    996              ENDIF
    997 
    998              r_a(j,i) = (pt_1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
     1000                pt1 = pt(k+1,j,i)
     1001                qv1 = q(k+1,j,i)
     1002             ENDIF
     1003
     1004             r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
     1005
     1006!
     1007!--          Make sure that the resistance does not drop to zero
     1008             IF ( ABS(r_a(j,i)) < 1.0E-10_wp )  r_a(j,i) = 1.0E-10_wp
    9991009
    10001010!
     
    10051015!--          night)
    10061016             IF ( radiation_scheme /= 'constant' )  THEN
    1007                 f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /&
    1008                               (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)        &
     1017                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /  &
     1018                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)          &
    10091019                               + 1.0_wp)) )
    10101020             ELSE
    10111021                f1 = 1.0_wp
    10121022             ENDIF
     1023
    10131024
    10141025!
     
    10401051!
    10411052!--             Calculate vapour pressure
    1042                 e  = qv_1 * surface_pressure / 0.622_wp
     1053                e  = qv1 * surface_pressure / 0.622_wp
    10431054                f3 = EXP ( -g_d(j,i) * (e_s - e) )
    10441055             ELSE
     
    10821093!--          In case of dewfall, set evapotranspiration to zero
    10831094!--          All super-saturated water is then removed from the air
    1084              IF ( humidity .AND. q_s <= qv_1 )  THEN
     1095             IF ( humidity .AND. q_s <= qv1 )  THEN
    10851096                r_canopy(j,i) = 0.0_wp
    10861097                r_soil(j,i)   = 0.0_wp
     
    11081119!--          reservoir are not allowed to take up water from the super-saturated
    11091120!--          air.
    1110              IF ( humidity )  THEN
    1111                 IF ( q_s <= qv_1 )  THEN
    1112                    IF ( .NOT. dewfall )  THEN
    1113                       f_qsws_veg  = 0.0_wp
    1114                       f_qsws_soil = 0.0_wp
    1115                       f_qsws_liq  = 0.0_wp
    1116                    ENDIF
    1117                 ENDIF
     1121             IF ( humidity .AND. q_s <= qv1 .AND. .NOT. dewfall )  THEN
     1122                f_qsws_veg  = 0.0_wp
     1123                f_qsws_soil = 0.0_wp
     1124                f_qsws_liq  = 0.0_wp
    11181125             ENDIF
    11191126
     
    11331140             rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
    11341141
    1135 
    11361142             IF ( humidity )  THEN
    11371143#if defined ( __rrtmg )
    11381144!
    11391145!--             Numerator of the prognostic equation
    1140                 coef_1 = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb)                &
     1146                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
    11411147                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
    1142                          + f_shf * pt_1 + f_qsws * ( qv_1 - q_s                &
     1148                         + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
    11431149                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
    11441150                         * t_soil(nzb_soil,j,i)
     
    11461152!
    11471153!--             Denominator of the prognostic equation
    1148                 coef_2 = rrtm_lwuflx_dt(0,nzb) + f_qsws * dq_s_dt              &
     1154                coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
    11491155                         + lambda_surface + f_shf / exn
    11501156#else
     
    11541160                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
    11551161                         * t_surface(j,i) ** 4                                 &
    1156                          + f_shf * pt_1 + f_qsws * ( qv_1                      &
     1162                         + f_shf * pt1 + f_qsws * ( qv1                        &
    11571163                         - q_s + dq_s_dt * t_surface(j,i) )                    &
    11581164                         + lambda_surface * t_soil(nzb_soil,j,i)
    11591165
    11601166!
    1161 !--                Denominator of the prognostic equation
    1162                    coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws   &
    1163                             * dq_s_dt + lambda_surface + f_shf / exn
     1167!--             Denominator of the prognostic equation
     1168                coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
     1169                         * dq_s_dt + lambda_surface + f_shf / exn
    11641170 
    11651171#endif
     
    11691175!
    11701176!--             Numerator of the prognostic equation
    1171                 coef_1 = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb)                &
     1177                coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
    11721178                         * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
    1173                          + f_shf * pt_1  + lambda_surface                      &
     1179                         + f_shf * pt1  + lambda_surface                       &
    11741180                         * t_soil(nzb_soil,j,i)
    11751181
    11761182!
    11771183!--             Denominator of the prognostic equation
    1178                 coef_2 = rrtm_lwuflx_dt(0,nzb) + lambda_surface + f_shf / exn
     1184                coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
    11791185#else
    11801186
     
    11821188!--             Numerator of the prognostic equation
    11831189                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
    1184                           * t_surface(j,i) ** 4 + f_shf * pt_1                 &
     1190                          * t_surface(j,i) ** 4 + f_shf * pt1                  &
    11851191                          + lambda_surface * t_soil(nzb_soil,j,i)
    11861192
     
    11911197#endif
    11921198             ENDIF
    1193 
    11941199
    11951200             tend = 0.0_wp
     
    12381243!--          Calculate fluxes
    12391244#if defined ( __rrtmg )
    1240              rad_net_l(j,i)   = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb)         &
     1245             rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
    12411246                                * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
    1242                                 - rrtm_lwuflx_dt(0,nzb) * t_surface_p(j,i)
     1247                                - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
    12431248
    12441249             IF ( rrtm_idrv == 1 )  THEN
    12451250                rad_net(j,i) = rad_net_l(j,i)
    12461251                rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
    1247                                       + rrtm_lwuflx_dt(0,nzb)                  &
     1252                                      + rad_lw_out_change_0(j,i)               &
    12481253                                      * ( t_surface_p(j,i) - t_surface(j,i) )
    12491254             ENDIF
     
    12521257                                * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
    12531258                                * t_surface(j,i)**3 * t_surface_p(j,i)
    1254              rad_net(j,i) = rad_net_l(j,i)
    12551259#endif
     1260
    12561261
    12571262             ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
    12581263                              - t_soil(nzb_soil,j,i))
    12591264
    1260              shf_eb(j,i)    = - f_shf * ( pt_1 - pt(k,j,i) )
     1265             shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
    12611266
    12621267             shf(j,i) = shf_eb(j,i) / rho_cp
    12631268
    12641269             IF ( humidity )  THEN
    1265                 qsws_eb(j,i)  = - f_qsws    * ( qv_1 - q_s + dq_s_dt           &
     1270                qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
    12661271                                * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
    12671272
    12681273                qsws(j,i) = qsws_eb(j,i) / rho_lv
    12691274
    1270                 qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv_1 - q_s               &
     1275                qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
    12711276                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
    12721277                                    * t_surface_p(j,i) )
    12731278
    1274                 qsws_soil_eb(j,i) = - f_qsws_soil * ( qv_1 - q_s               &
     1279                qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
    12751280                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
    12761281                                    * t_surface_p(j,i) )
    12771282
    1278                 qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv_1 - q_s               &
     1283                qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
    12791284                                    + dq_s_dt * t_surface(j,i) - dq_s_dt       &
    12801285                                    * t_surface_p(j,i) )
     
    12851290                r_s(j,i) = 1.0E10_wp
    12861291             ELSE
    1287                 r_s(j,i) = - rho_lv * ( qv_1 - q_s + dq_s_dt                   &
     1292                r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                   &
    12881293                           * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
    12891294                           / qsws_eb(j,i) - r_a(j,i)
  • palm/trunk/SOURCE/radiation_model.f90

    r1702 r1709  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
     22! corrections
    2223!
    2324! Former revisions:
     
    197198    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
    198199                alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
     200                rad_lw_out_change_0,         & !< change in LW out due to change in surface temperature
    199201                rad_net,                     & !< net radiation at the surface
    200202                rad_net_av                     !< average of rad_net
     
    363365           radiation_clearsky, radiation_rrtmg, radiation_scheme,              &
    364366           radiation_tendency, rad_lw_in, rad_lw_in_av, rad_lw_out,            &
    365            rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,            &
    366            rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av,   &
    367            rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,   &
    368            skip_time_do_radiation, sw_radiation, time_radiation, time_utc_init
     367           rad_lw_out_av, rad_lw_out_change_0, rad_lw_cs_hr, rad_lw_cs_hr_av,  &
     368           rad_lw_hr, rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
     369           rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,            &
     370           rad_sw_hr_av, sigma_sb, skip_time_do_radiation, sw_radiation,       &
     371           time_radiation, time_utc_init
    369372
    370373
    371374#if defined ( __rrtmg )
    372     PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv,          &
    373            rrtm_lwuflx_dt
     375    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv
    374376#endif
    375377
     
    390392          ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
    391393          rad_net = 0.0_wp
     394       ENDIF
     395
     396!
     397!--    Allocate array for storing the surface net radiation
     398       IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) )  THEN
     399          ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
     400          rad_lw_out_change_0 = 0.0_wp
    392401       ENDIF
    393402
     
    654663       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    655664          IF ( radiation_scheme == 'clear-sky' )  THEN
    656               CALL radiation_clearsky
     665             CALL radiation_clearsky
    657666          ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
    658667             CALL radiation_rrtmg
     
    679688       INTEGER(iwp) :: i, j, k   !< loop indices
    680689       REAL(wp)     :: exn,   &  !< Exner functions at surface
    681                        exn_1, &  !< Exner functions at first grid level
    682                        pt_1      !< potential temperature at first grid level
     690                       exn1, &  !< Exner functions at first grid level
     691                       pt1       !< potential temperature at first grid level
    683692
    684693!
     
    696705!--    Calculate radiation fluxes and net radiation (rad_net) for each grid
    697706!--    point
    698        DO i = nxl, nxr
    699           DO j = nys, nyn
     707       DO i = nxlg, nxrg
     708          DO j = nysg, nyng
    700709             k = nzb_s_inner(j,i)
    701710
    702              exn_1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
     711             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
    703712
    704713             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
     
    707716
    708717             IF ( cloud_physics )  THEN
    709                 pt_1 = pt(k+1,j,i) + l_d_cp / exn_1 * ql(k+1,j,i)
    710                 rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt_1 * exn_1)**4
     718                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
     719                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
    711720             ELSE
    712                 rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn_1)**4
     721                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
    713722             ENDIF
    714723
     
    718727          ENDDO
    719728       ENDDO
    720 
    721        CALL exchange_horiz_2d( rad_lw_in,  nbgp )
    722        CALL exchange_horiz_2d( rad_lw_out, nbgp )
    723        CALL exchange_horiz_2d( rad_sw_in,  nbgp )
    724        CALL exchange_horiz_2d( rad_sw_out, nbgp )
    725        CALL exchange_horiz_2d( rad_net,    nbgp )
    726 
    727        RETURN
    728729
    729730    END SUBROUTINE radiation_clearsky
     
    901902                   rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k) * d_hours_day
    902903                ENDDO
     904
     905!
     906!--             Save change in LW heating rate
     907                rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
    903908
    904909             ENDIF
     
    953958
    954959       CALL exchange_horiz_2d( rad_net, nbgp )
     960       CALL exchange_horiz_2d( rad_lw_out_change_0, nbgp )
    955961#endif
    956962
     
    13471353       ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
    13481354
    1349        rrtm_lwuflxc_dt = 0.0_wp
     1355       rrtm_lwuflx_dt = 0.0_wp
    13501356       rrtm_lwuflxc_dt = 0.0_wp
    13511357
  • palm/trunk/SOURCE/read_3d_binary.f90

    r1692 r1709  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Added rad_lw_out_change_0, increased binary_version
    2222!
    2323! Former revisions:
     
    124124        ONLY: rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,        &
    125125              rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av,          &
    126               rad_lw_out_av, rad_sw_in, rad_sw_in_av, rad_sw_out,              &
    127               rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,         &
    128               rad_sw_hr_av
     126              rad_lw_out_av, rad_lw_out_change_0, rad_sw_in, rad_sw_in_av,     &
     127              rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av,        &
     128              rad_sw_hr, rad_sw_hr_av
    129129
    130130    USE random_function_mod,                                                   &
     
    348348!--    First compare the version numbers
    349349       READ ( 13 )  version_on_file
    350        binary_version = '4.1'
     350       binary_version = '4.2'
    351351       IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
    352352          WRITE( message_string, * ) 'version mismatch concerning data ',      &
     
    877877                   rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    878878                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     879
     880                CASE ( 'rad_lw_out_change_0' )
     881                   IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
     882                      ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
     883                   ENDIF
     884                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     885                   rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
     886                              = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    879887
    880888                CASE ( 'rad_lw_cs_hr' )
  • palm/trunk/SOURCE/surface_layer_fluxes.f90

    r1706 r1709  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Bugfix: division by zero could occur when calculating rib at low wind speeds
     22! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
     23! neutral = .T.
    2224!
    2325! Former revisions:
     
    9698!>    that this method cannot be used in case of roughness heterogeneity
    9799!>
    98 !> @todo limiting of uv_total might be required in some cases
    99100!> @todo (re)move large_scale_forcing actions
    100101!> @todo check/optimize OpenMP and OpenACC directives
     
    208209          CALL calc_scaling_parameters
    209210
     211          CALL calc_uv_total
     212
    210213          IF ( .NOT. neutral )  THEN
    211214             CALL calc_ol
     
    220223!--    number to calculate the Obukhov length
    221224       ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' )  THEN
     225
     226          CALL calc_uv_total
    222227
    223228          IF ( .NOT. neutral )  THEN
     
    273278!--    Allocate field for storing the horizontal velocity
    274279       ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) )
     280
     281
     282!
     283!--    In case of runs with neutral statification, set Obukhov length to a
     284!--    large value
     285       IF ( neutral ) ol = 1.0E10_wp
    275286
    276287       IF ( most_method == 'lookup' )  THEN
     
    384395! Description:
    385396! ------------
     397!> Compute the absolute value of the horizontal velocity (relative to the
     398!> surface). This is required by all methods
     399!------------------------------------------------------------------------------!
     400    SUBROUTINE calc_uv_total
     401
     402       IMPLICIT NONE
     403
     404
     405       !$OMP PARALLEL DO PRIVATE( k, z_mo )
     406       !$acc kernels loop
     407       DO  i = nxl, nxr
     408          DO  j = nys, nyn
     409
     410             k   = nzb_s_inner(j,i)
     411             uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
     412                                         - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
     413                              ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
     414                                         - v(k,j,i)   - v(k,j+1,i) ) )**2 )
     415
     416!
     417!--          For too small values of the local wind, MOST does not work. A
     418!--          threshold value is thus set if required
     419!            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
     420
     421          ENDDO
     422       ENDDO
     423
     424!
     425!--    Values of uv_total need to be exchanged at the ghost boundaries
     426       !$acc update host( uv_total )
     427       CALL exchange_horiz_2d( uv_total )
     428       !$acc update device( uv_total )
     429
     430    END SUBROUTINE calc_uv_total
     431
     432
     433!------------------------------------------------------------------------------!
     434! Description:
     435! ------------
    386436!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
    387437!------------------------------------------------------------------------------!
     
    394444
    395445       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number
    396                        
     446
    397447       REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
    398448                       f_d_ol, & !< Derivative of f
    399449                       ol_l,   & !< Lower bound of L for Newton iteration
    400450                       ol_m,   & !< Previous value of L for Newton iteration
    401                        ol_old, & !< Previous time step value of L 
     451                       ol_old, & !< Previous time step value of L
    402452                       ol_u      !< Upper bound of L for Newton iteration
    403453
    404 !
    405 !--    Compute the absolute value of the horizontal velocity (relative to the
    406 !--    surface). This is required by all methods
    407        !$OMP PARALLEL DO PRIVATE( k, z_mo )
    408        !$acc kernels loop
    409        DO  i = nxl, nxr
    410           DO  j = nys, nyn
    411 
    412              k   = nzb_s_inner(j,i)
    413 
    414              uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)      &
    415                                          - u(k,j,i)   - u(k,j,i+1) ) )**2 +    &
    416                               ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)           &
    417                                          - v(k,j,i)   - v(k,j+1,i) ) )**2 )
    418 
    419 !
    420 !--          For too small values of the local wind, MOST does not work. A
    421 !--          threshold value is thus set if required
    422 !            uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
    423 
    424           ENDDO
    425        ENDDO
    426 
    427 !
    428 !--    Values of uv_total need to be exchanged at the ghost boundaries
    429        !$acc update host( uv_total )
    430        CALL exchange_horiz_2d( uv_total )
    431        !$acc update device( uv_total )
    432454
    433455       IF ( TRIM( most_method ) /= 'circular' )  THEN
     
    447469                   IF ( humidity )  THEN
    448470                      rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) )      &
    449                            / ( uv_total(j,i)**2 * vpt(k+1,j,i) )
     471                           / ( uv_total(j,i)**2 * vpt(k+1,j,i) + 1.0E-20_wp )
    450472                   ELSE
    451473                      rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) )        &
    452                            / ( uv_total(j,i)**2 * pt(k+1,j,i) )
     474                           / ( uv_total(j,i)**2 * pt(k+1,j,i)  + 1.0E-20_wp )
    453475                   ENDIF     
    454476                ELSE
     
    462484                                 * q(k+1,j,i) ) * shf(j,i) + 0.61_wp           &
    463485                                 * pt(k+1,j,i) * qsws(j,i) )   &
    464                                  / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2 )
     486                                 / ( uv_total(j,i)**3 * vpt(k+1,j,i) * kappa**2&
     487                                 + 1.0E-20_wp)
    465488                   ELSE
    466489                      rib(j,i) = - g * z_mo * shf(j,i)                         &
    467                            / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2 )
     490                           / ( uv_total(j,i)**3 * pt(k+1,j,i) * kappa**2       &
     491                           + 1.0E-20_wp )
    468492                   ENDIF
    469493                ENDIF 
  • palm/trunk/SOURCE/write_3d_binary.f90

    r1692 r1709  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Added rad_lw_out_change_0, increased binary_version
    2222!
    2323! Former revisions:
     
    114114    USE radiation_model_mod,                                                   &
    115115        ONLY: radiation, rad_net, rad_net_av, rad_lw_in, rad_lw_in_av,         &
    116               rad_lw_out, rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av,        &
    117               rad_lw_hr, rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out,    &
    118               rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,         &
    119               rad_sw_hr_av
     116              rad_lw_out, rad_lw_out_av, rad_lw_out_change_0, rad_lw_cs_hr,   &
     117              rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,             &
     118              rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,           &
     119              rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av
    120120
    121121    USE random_function_mod,                                                   &
     
    139139!
    140140!-- Write arrays.
    141     binary_version = '4.1'
     141    binary_version = '4.2'
    142142
    143143    WRITE ( 14 )  binary_version
     
    297297          WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av 
    298298       ENDIF
     299       IF ( ALLOCATED( rad_lw_out_change_0 ) )  THEN
     300          WRITE ( 14 )  'rad_lw_out_change_0 '
     301          WRITE ( 14 )  rad_lw_out_change_0
     302       ENDIF
    299303       IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
    300304          WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
Note: See TracChangeset for help on using the changeset viewer.