Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Jun 12, 2019 11:52:39 AM (5 years ago)
Author:
suehring
Message:

Synthetic turbulence generator: Revise bias correction of imposed perturbations (correction via volume flow can create instabilities in case the mean volume flow is close to zero); Introduce lower limits in calculation of coefficient matrix, else the calculation may become numerically unstable; Impose perturbations every timestep, even though no new set of perturbations is generated in case dt_stg_call /= dt_3d; Implement a gradual decrease of Reynolds stress and length scales above ABL height (within 1 length scale above ABL depth to 1/10) rather than an abrupt decline; Bugfix in non-nested case: use ABL height for parametrized turbulence; Offline nesting: Rename subroutine for ABL height calculation and make it public; Change top boundary condition for pressure back to Neumann

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/nesting_offl_mod.f90

    r3987 r4022  
    2525! -----------------
    2626! $Id$
     27! Detection of boundary-layer depth in stable boundary layer on basis of
     28! boundary data improved
     29! Routine for boundary-layer depth calculation renamed and made public
     30!
     31! 3987 2019-05-22 09:52:13Z kanani
    2732! Introduce alternative switch for debug output during timestepping
    2833!
    … …  
    117122!
    118123!-- Public subroutines
    119     PUBLIC nesting_offl_bc, nesting_offl_check_parameters, nesting_offl_header,&
    120            nesting_offl_init, nesting_offl_mass_conservation, nesting_offl_parin
     124    PUBLIC nesting_offl_bc,                                                    &
     125           nesting_offl_calc_zi,                                               &
     126           nesting_offl_check_parameters,                                      &
     127           nesting_offl_header,                                                &
     128           nesting_offl_init,                                                  &
     129           nesting_offl_mass_conservation,                                     &
     130           nesting_offl_parin
    121131!
    122132!-- Public variables
    … …  
    126136       MODULE PROCEDURE nesting_offl_bc
    127137    END INTERFACE nesting_offl_bc
     138   
     139    INTERFACE nesting_offl_calc_zi
     140       MODULE PROCEDURE nesting_offl_calc_zi
     141    END INTERFACE nesting_offl_calc_zi
    128142   
    129143    INTERFACE nesting_offl_check_parameters
    … …  
    806820!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
    807821!--    Therefore, calculate boundary-layer depth first.
    808        CALL calc_zi
     822       CALL nesting_offl_calc_zi
    809823       CALL adjust_sponge_layer
    810824   
    … …  
    833847!> bulk-Richardson criterion.
    834848!------------------------------------------------------------------------------!
    835     SUBROUTINE calc_zi
     849    SUBROUTINE nesting_offl_calc_zi
    836850       
    837851       USE basic_constants_and_equations_mod,                                  &
    … …  
    845859       IMPLICIT NONE
    846860
    847        INTEGER(iwp) :: i            !< loop index in x-direction
    848        INTEGER(iwp) :: j            !< loop index in y-direction
    849        INTEGER(iwp) :: k            !< loop index in z-direction
    850        INTEGER(iwp) :: k_surface    !< topography top index in z-direction
     861       INTEGER(iwp) :: i                            !< loop index in x-direction
     862       INTEGER(iwp) :: j                            !< loop index in y-direction
     863       INTEGER(iwp) :: k                            !< loop index in z-direction
     864       INTEGER(iwp) :: k_max_loc                    !< index of maximum wind speed along z-direction
     865       INTEGER(iwp) :: k_surface                    !< topography top index in z-direction
     866       INTEGER(iwp) :: num_boundary_gp_non_cyclic   !< number of non-cyclic boundaries, used for averaging ABL depth
     867       INTEGER(iwp) :: num_boundary_gp_non_cyclic_l !< number of non-cyclic boundaries, used for averaging ABL depth
    851868   
    852869       REAL(wp) ::  ri_bulk                 !< bulk Richardson number
    … …  
    854871       REAL(wp) ::  topo_max                !< maximum topography level in model domain
    855872       REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
    856        REAL(wp) ::  u_comp                  !< u-component
    857        REAL(wp) ::  v_comp                  !< v-component
    858873       REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
    859874       REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
    … …  
    861876   
    862877       REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
     878       REAL(wp), DIMENSION(nzb:nzt+1) ::  uv_abs  !< vertical profile of horizontal wind speed at (j,i)-grid point
    863879
    864880       
    … …  
    867883!--    Start with the left and right boundaries.
    868884       zi_l      = 0.0_wp
     885       num_boundary_gp_non_cyclic_l = 0
    869886       IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
     887!
     888!--       Sum-up and store number of boundary grid points used for averaging
     889!--       ABL depth
     890          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
     891                                         nxr - nxl + 1
    870892!
    871893!--       Determine index along x. Please note, index indicates boundary
    … …  
    898920!--          the boundary, where velocity inside the building is zero
    899921!--          (k_surface is the index of the lowest upward-facing surface).
     922             uv_abs(:) = SQRT( MERGE( u(:,j,i+1), u(:,j,i),                   &
     923                                      bc_dirichlet_l )**2 +                   &
     924                               v(:,j,i)**2 )
     925!
     926!--          Determine index of the maximum wind speed                               
     927             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
     928
    900929             zi_local = 0.0_wp
    901930             DO  k = k_surface+1, nzt
    902                 u_comp = MERGE( u(k,j,i+1), u(k,j,i), bc_dirichlet_l )
    903                 v_comp = v(k,j,i)
    904931                ri_bulk = zu(k) * g / vpt_surface *                            &
    905932                          ( vpt_col(k) - vpt_surface ) /                       &
    906                           ( u_comp * u_comp + v_comp * v_comp + 1E-5_wp )
    907                        
    908                 IF ( zi_local == 0.0_wp  .AND.  ri_bulk > ri_bulk_crit )       &
     933                          ( uv_abs(k) + 1E-5_wp )
     934!
     935!--             Check if critical Richardson number is exceeded. Further, check
     936!--             if there is a maxium in the wind profile in order to detect also
     937!--             ABL heights in the stable boundary layer.
     938                IF ( zi_local == 0.0_wp  .AND.                                 &
     939                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
    909940                   zi_local = zu(k)           
    910941             ENDDO
    … …  
    920951!--    Do the same at the north and south boundaries.
    921952       IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
     953       
     954          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
     955                                         nxr - nxl + 1
    922956       
    923957          j = MERGE( -1, nyn + 1, bc_dirichlet_s )
    … …  
    935969             ENDIF
    936970         
     971             uv_abs(:) = SQRT( u(:,j,i)**2 +                                   &
     972                               MERGE( v(:,j+1,i), v(:,j,i),                    &
     973                               bc_dirichlet_s )**2 )
     974!
     975!--          Determine index of the maximum wind speed
     976             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
     977         
    937978             zi_local = 0.0_wp
    938979             DO  k = k_surface+1, nzt               
    939                 u_comp = u(k,j,i)
    940                 v_comp = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s )
    941980                ri_bulk = zu(k) * g / vpt_surface *                            &
    942981                       ( vpt_col(k) - vpt_surface ) /                          &
    943                        ( u_comp * u_comp + v_comp * v_comp + 1E-5_wp )
    944                        
    945                 IF ( zi_local == 0.0_wp  .AND.  ri_bulk > 0.25_wp )            &
     982                       ( uv_abs(k) + 1E-5_wp )
     983!
     984!--             Check if critical Richardson number is exceeded. Further, check
     985!--             if there is a maxium in the wind profile in order to detect also
     986!--             ABL heights in the stable boundary layer.                       
     987                IF ( zi_local == 0.0_wp  .AND.                                 &
     988                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
    946989                   zi_local = zu(k)   
    947990             ENDDO
    … …  
    955998       CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
    956999                           comm2d, ierr )
     1000       CALL MPI_ALLREDUCE( num_boundary_gp_non_cyclic_l,                       &
     1001                           num_boundary_gp_non_cyclic,                         &
     1002                           1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
    9571003#else
    9581004       zi_ribulk = zi_l
     1005       num_boundary_gp_non_cyclic = num_boundary_gp_non_cyclic_l
    9591006#endif
    960        zi_ribulk = zi_ribulk / REAL( 2 * nx + 2 * ny, KIND = wp )
     1007       zi_ribulk = zi_ribulk / REAL( num_boundary_gp_non_cyclic, KIND = wp )
    9611008!
    9621009!--    Finally, check if boundary layer depth is not below the any topography.
    … …  
    9681015       
    9691016#if defined( __parallel )
    970        CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,            &
     1017       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,         &
    9711018                           comm2d, ierr )
    9721019#else
    9731020       topo_max     = topo_max_l
    9741021#endif
    975        zi_ribulk = MAX( zi_ribulk, topo_max )
    976        
    977     END SUBROUTINE calc_zi
     1022!        zi_ribulk = MAX( zi_ribulk, topo_max )
     1023       
     1024    END SUBROUTINE nesting_offl_calc_zi
    9781025   
    9791026   
    … …  
    13191366!--    boundary data. This is requiered for initialize the synthetic turbulence
    13201367!--    generator correctly.
    1321        CALL calc_zi
     1368       CALL nesting_offl_calc_zi
    13221369       
    13231370!
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |