Ignore:
Timestamp:
Jun 19, 2018 2:03:12 PM (3 years ago)
Author:
gronemeier
Message:

merge with branch rans

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/model_1d_mod.f90

    r3049 r3083  
    2525! -----------------
    2626! $Id$
     27! Bugfixes:
     28!   - preset te_diss and te_e to avoid runtime errors
     29!   - implementation of buoyancy term to dissipation
     30!     according to Sogachev et al. (2012)
     31!   - where diss_p < 0 set diss_p = 0.1 diss
     32!   - calculate progn eq(diss) starting from nzb+1
     33! Changes:
     34!   - add sig_e to TKE equation
     35!   - adjust prognostic equation of diss
     36!   - set model constants according to Koblitz (2013)
     37!   - renamed c_m to c_0
     38!   - rename l_black into l1d_init
     39!   - calculate l_grid within init_1d_model and save it as l1d_init
     40!   - calculate l1d according to DE85 if dissipation is a prognostic value
     41!   - made annotations doxygen-readable
     42!
     43! 3049 2018-05-29 13:52:36Z Giersch
    2744! Error messages revised
    2845!
     
    164181    USE kinds
    165182
    166     USE pegrid
     183    USE pegrid,                                                                &
     184        ONLY:  myid
    167185       
    168186
     
    175193    LOGICAL ::  stop_dt_1d = .FALSE.             !< termination flag, used in case of too small timestep (1d-model)
    176194
    177     REAL(wp) ::  c_1 = 1.44_wp                 !< model constant
    178     REAL(wp) ::  c_2 = 1.92_wp                 !< model constant
    179     REAL(wp) ::  c_3 = 1.44_wp                 !< model constant
    180     REAL(wp) ::  c_h = 0.0015_wp               !< model constant according to Detering and Etling (1985)
    181     REAL(wp) ::  c_m = 0.4_wp                  !< model constant, 0.4 according to Detering and Etling (1985)
    182     REAL(wp) ::  c_mu = 0.09_wp                !< model constant
     195    REAL(wp) ::  alpha_buoyancy                !< model constant according to Koblitz (2013)
     196    REAL(wp) ::  c_0 = 0.03_wp**0.25_wp        !< model constant according to Koblitz (2013)
     197    REAL(wp) ::  c_1 = 1.52_wp                 !< model constant according to Koblitz (2013)
     198    REAL(wp) ::  c_2 = 1.83_wp                 !< model constant according to Koblitz (2013)
     199    REAL(wp) ::  c_3                           !< model constant
     200    REAL(wp) ::  c_mu                          !< model constant
    183201    REAL(wp) ::  damp_level_1d = -1.0_wp       !< namelist parameter
    184202    REAL(wp) ::  dt_1d = 60.0_wp               !< dynamic timestep (1d-model)
     
    187205    REAL(wp) ::  dt_run_control_1d = 60.0_wp   !< namelist parameter
    188206    REAL(wp) ::  end_time_1d = 864000.0_wp     !< namelist parameter
     207    REAL(wp) ::  lambda                        !< maximum mixing length
    189208    REAL(wp) ::  qs1d                          !< characteristic humidity scale (1d-model)
    190209    REAL(wp) ::  simulated_time_1d = 0.0_wp    !< updated simulated time (1d-model)
    191     REAL(wp) ::  sig_diss = 1.3_wp             !< model constant
     210    REAL(wp) ::  sig_diss = 2.95_wp            !< model constant according to Koblitz (2013)
     211    REAL(wp) ::  sig_e = 2.95_wp               !< model constant according to Koblitz (2013)
    192212    REAL(wp) ::  time_pr_1d = 0.0_wp           !< updated simulated time for profile output (1d-model)
    193213    REAL(wp) ::  time_run_control_1d = 0.0_wp  !< updated simulated time for run-control output (1d-model)
     
    198218    REAL(wp) ::  z01d                          !< roughness length for momentum (1d-model)
    199219    REAL(wp) ::  z0h1d                         !< roughness length for scalars (1d-model)
    200 
    201220
    202221    REAL(wp), DIMENSION(:), ALLOCATABLE ::  diss1d   !< tke dissipation rate (1d-model)
     
    279298
    280299    INTEGER(iwp) ::  k  !< loop index
    281    
    282     REAL(wp) ::  lambda !< maximum mixing length
    283300
    284301!
     
    324341!
    325342!--       Use the same mixing length as in 3D model (LES-mode)
    326           !@todo: rename (delete?) this option
    327           ! As the mixing length is different between RANS and LES mode, it
    328           ! must be distinguished here between these modes. For RANS mode,
    329           ! the mixing length is calculated accoding to Blackadar, which is
    330           ! the other option at this point.
    331           ! Maybe delete this option entirely (not appropriate in LES case)
    332           ! 2018-03-20, gronemeier
     343          !> @todo rename (delete?) this option
     344          !> As the mixing length is different between RANS and LES mode, it
     345          !> must be distinguished here between these modes. For RANS mode,
     346          !> the mixing length is calculated accoding to Blackadar, which is
     347          !> the other option at this point.
     348          !> Maybe delete this option entirely (not appropriate in LES case)
     349          !> 2018-03-20, gronemeier
    333350          DO  k = nzb+1, nzt
    334351             l1d_init(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
     
    359376       us1d = 0.1_wp   ! without initial friction the flow would not change
    360377    ELSE
    361        diss1d(nzb+1) = 1.0_wp
     378       diss1d(nzb+1) = 0.001_wp
    362379       e1d(nzb+1)  = 1.0_wp
    363380       km1d(nzb+1) = 1.0_wp
     
    372389
    373390!
    374 !-- Tendencies must be preset in order to avoid runtime errors within the
    375 !-- first Runge-Kutta step
     391!-- Tendencies must be preset in order to avoid runtime errors
     392    te_diss  = 0.0_wp
    376393    te_dissm = 0.0_wp
     394    te_e  = 0.0_wp
    377395    te_em = 0.0_wp
    378396    te_um = 0.0_wp
     
    381399!
    382400!-- Set model constant
    383     IF ( dissipation_1d == 'as_in_3d_model' )  c_m = 0.1_wp
     401    IF ( dissipation_1d == 'as_in_3d_model' )  c_0 = 0.1_wp
     402    c_mu = c_0**4
    384403
    385404!
     
    423442!-- Determine the time step at the start of a 1D-simulation and
    424443!-- determine and printout quantities used for run control
    425     dt_1d = 1.0_wp
     444    dt_1d = 0.01_wp
    426445    CALL run_control_1d
    427446
     
    483502!--             dissipation rate
    484503                IF ( dissipation_1d == 'detering' )  THEN
    485                    diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
     504                   diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    486505                ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    487506                   diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) &
     
    497516                                     kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
    498517                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
    499                                                  ) * ddzw(k)                   &
     518                                                 ) * ddzw(k) / sig_e           &
    500519                                   - diss1d(k)
    501520
     
    503522!
    504523!--                dissipation rate
    505                    te_diss(k) = km1d(k) *                                      &
     524                   IF ( rif1d(k) >= 0.0_wp )  THEN
     525                      alpha_buoyancy = 1.0_wp - l1d(k) / lambda
     526                   ELSE
     527                      alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp )    &
     528                                                         / ( c_2 - c_1    ) )  &
     529                                              * l1d(k) / lambda
     530                   ENDIF
     531                   c_3 = ( c_1 - c_2 ) * alpha_buoyancy
     532                   te_diss(k) = ( km1d(k) *                                    &
    506533                                  ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2  &
    507534                                  + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2  &
    508                                   ) * c_1 * c_mu**0.75 / c_h * f / us1d        &
    509                                     * SQRT(e1d(k))                             &
     535                                  ) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda )  &
    510536                                  - g / pt_0 * kh1d(k) * flux * c_3            &
    511                                     * diss1d(k) / ( e1d(k) + 1.0E-20_wp )      &
    512                                   + ( kmzp * ( diss1d(k+1) - diss1d(k) )       &
     537                                  - c_2 * diss1d(k)                            &
     538                                ) * diss1d(k) / ( e1d(k) + 1.0E-20_wp )        &
     539                                + (   kmzp * ( diss1d(k+1) - diss1d(k) )       &
    513540                                           * ddzu(k+1)                         &
    514541                                    - kmzm * ( diss1d(k) - diss1d(k-1) )       &
    515542                                           * ddzu(k)                           &
    516                                     ) * ddzw(k) / sig_diss                     &
    517                                   - c_2 * diss1d(k)**2 / ( e1d(k) + 1.0E-20_wp )
     543                                  ) * ddzw(k) / sig_diss
    518544
    519545                ENDIF
     
    546572!--          dissipation rate
    547573             IF ( dissipation_1d == 'detering' )  THEN
    548                 diss1d(k) = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
     574                diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)
    549575             ELSEIF ( dissipation_1d == 'as_in_3d_model' )  THEN
    550576                diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) )  &
     
    565591!--          TKE
    566592             IF ( .NOT. dissipation_1d == 'prognostic' )  THEN
     593                !> @query why integrate over 2dz
     594                !>   Why is it allowed to integrate over two delta-z for e
     595                !>   while for u and v it is not?
     596                !>   2018-04-23, gronemeier
    567597                te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2&
    568598                                    + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2&
     
    572602                                     kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1)  &
    573603                                   - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k)    &
    574                                                  ) * ddzw(k)                   &
     604                                                 ) * ddzw(k) / sig_e           &
    575605                                   - diss1d(k)
    576606             ENDIF
     
    595625             ENDDO
    596626
    597              IF ( dissipation_1d == 'prognostic' )  THEN
    598                 DO  k = nzb_diff, nzt
    599                    diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
    600                                                  tsc(3) * te_dissm(k) )
    601                 ENDDO
    602              ENDIF
    603627!
    604628!--          Eliminate negative TKE values, which can result from the
     
    606630!--          value is reduced to 10 percent of its old value.
    607631             WHERE ( e1d_p < 0.0_wp )  e1d_p = 0.1_wp * e1d
     632
     633             IF ( dissipation_1d == 'prognostic' )  THEN
     634                DO  k = nzb+1, nzt
     635                   diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &
     636                                                       tsc(3) * te_dissm(k) )
     637                ENDDO
     638                WHERE ( diss1d_p < 0.0_wp )  diss1d_p = 0.1_wp * diss1d
     639             ENDIF
    608640          ENDIF
    609641
     
    643675                   IF ( dissipation_1d == 'prognostic' )  THEN
    644676                      DO k = nzb+1, nzt
    645                          te_dissm(k) = -9.5625_wp * te_diss(k) + 5.3125_wp * te_dissm(k)
     677                         te_dissm(k) = -9.5625_wp * te_diss(k)  &
     678                                     +  5.3125_wp * te_dissm(k)
    646679                      ENDDO
    647680                   ENDIF
     
    650683             ENDIF
    651684          ENDIF
    652 
    653685
    654686!
     
    701733
    702734             ENDIF    ! constant_flux_layer
    703 
     735             !> @todo combine if clauses
     736             !>   The previous and following if clauses can be combined into a
     737             !>   single clause
     738             !>   2018-04-23, gronemeier
    704739!
    705740!--          Compute the Richardson-flux numbers,
     
    789824!--             compatibility with the 3D model.
    790825                IF ( ibc_e_b == 2 )  THEN
    791                    e1d(nzb+1) = ( us1d / c_m )**2
     826                   e1d(nzb+1) = ( us1d / c_0 )**2
    792827                ENDIF
    793828                IF ( dissipation_1d == 'prognostic' )  THEN
    794                    e1d(nzb+1) = us1d**2 / SQRT( c_mu )
     829                   e1d(nzb+1) = ( us1d / c_0 )**2
    795830                   diss1d(nzb+1) = us1d**3 / ( kappa * zu(nzb+1) )
    796831                   diss1d(nzb) = diss1d(nzb+1)
     
    829864!--          in the dissipation of TKE via l1d_diss. Otherwise, km1d would be
    830865!--          too large.
    831              IF ( mixing_length_1d == 'blackadar' )  THEN
     866             IF ( dissipation_1d /= 'prognostic' )  THEN
     867                IF ( mixing_length_1d == 'blackadar' )  THEN
     868                   DO  k = nzb+1, nzt
     869                      IF ( rif1d(k) >= 0.0_wp )  THEN
     870                         l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
     871                         l1d_diss(k) = l1d(k)
     872                      ELSE
     873                         l1d(k) = l1d_init(k)
     874                         l1d_diss(k) = l1d_init(k) *                           &
     875                                       SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
     876                      ENDIF
     877                   ENDDO
     878                ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
     879                   DO  k = nzb+1, nzt
     880                      dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
     881                      IF ( dpt_dz > 0.0_wp )  THEN
     882                         l_stable = 0.76_wp * SQRT( e1d(k) )                   &
     883                                  / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
     884                      ELSE
     885                         l_stable = l1d_init(k)
     886                      ENDIF
     887                      l1d(k) = MIN( l1d_init(k), l_stable )
     888                      l1d_diss(k) = l1d(k)
     889                   ENDDO
     890                ENDIF
     891             ELSE
    832892                DO  k = nzb+1, nzt
    833                    IF ( rif1d(k) >= 0.0_wp )  THEN
    834                       l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) )
    835                       l1d_diss(k) = l1d(k)
    836                    ELSE
    837                       l1d(k) = l1d_init(k)
    838                       l1d_diss(k) = l1d_init(k) *                              &
    839                                     SQRT( 1.0_wp - 16.0_wp * rif1d(k) )
    840                    ENDIF
    841                 ENDDO
    842              ELSEIF ( mixing_length_1d == 'as_in_3d_model' )  THEN
    843                 DO  k = nzb+1, nzt
    844                    dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k)
    845                    IF ( dpt_dz > 0.0_wp )  THEN
    846                       l_stable = 0.76_wp * SQRT( e1d(k) ) /                    &
    847                                      SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp
    848                    ELSE
    849                       l_stable = l1d_init(k)
    850                    ENDIF
    851                    l1d(k) = MIN( l1d_init(k), l_stable )
    852                    l1d_diss(k) = l1d(k)
     893                   l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) )                   &
     894                          / ( diss1d(k) + 1.0E-30_wp )
    853895                ENDDO
    854896             ENDIF
     
    867909                ENDIF
    868910             ENDIF
     911
    869912             IF ( dissipation_1d == 'prognostic' )  THEN
    870913                DO  k = nzb_diff, nzt
     
    873916             ELSE
    874917                DO  k = nzb_diff, nzt
    875                    km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k)
     918                   km1d(k) = c_0 * SQRT( e1d(k) ) * l1d(k)
    876919                ENDDO
    877920             ENDIF
     
    10181061
    10191062    REAL(wp) ::  dt_diff  !< time step accorind to diffusion criterion
     1063    REAL(wp) ::  dt_old   !< previous time step
    10201064    REAL(wp) ::  fac      !< factor of criterion
    10211065    REAL(wp) ::  value    !< auxiliary variable
    10221066
    10231067!
     1068!-- Save previous time step
     1069    dt_old = dt_1d
     1070
     1071!
    10241072!-- Compute the currently feasible time step according to the diffusion
    10251073!-- criterion. At nzb+1 the half grid length is used.
    1026     fac = 0.125  !0.35_wp                                                       !### changed from 0.35
     1074    fac = 0.125
    10271075    dt_diff = dt_max_1d
    10281076    DO  k = nzb+2, nzt
     
    10341082
    10351083!
     1084!-- Limit the new time step to a maximum of 10 times the previous time step
     1085    dt_1d = MIN( dt_old * 10.0_wp, dt_1d )
     1086
     1087!
    10361088!-- Set flag when the time step becomes too small
    1037     IF ( dt_1d < ( 0.00001_wp * dt_max_1d ) )  THEN
     1089    IF ( dt_1d < ( 1.0E-15_wp * dt_max_1d ) )  THEN
    10381090       stop_dt_1d = .TRUE.
    10391091
Note: See TracChangeset for help on using the changeset viewer.