Changeset 4340


Ignore:
Timestamp:
Dec 16, 2019 8:17:03 AM (4 years ago)
Author:
Giersch
Message:

Topography closed channel flow with symmetric boundaries implemented, ID tag in radiation module corrected

Location:
palm/trunk/SOURCE
Files:
5 edited

Legend:

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

    r4330 r4340  
    2525! -----------------
    2626! $Id$
     27! Topography closed channel flow with symmetric boundaries implemented
     28!
     29! 4330 2019-12-10 16:16:33Z knoop
    2730! Bugix: removed syntax error introduced by last commit
    2831!
     
    166169               salsa,                                                          &
    167170               scalar_advec,                                                   &
     171               symmetry_flag,                                                  &
    168172               intermediate_timestep_count,                                    &
    169173               u_gtrans,                                                       &
     
    502506                             BTEST(wall_flags_static_0(k+1,j,i),1)       .AND.     &
    503507                             BTEST(wall_flags_static_0(k,j,i),1) )       .OR.      &
    504                      k == nzt )                                                    &
     508                     ( k == nzt .AND. symmetry_flag == 0 ) )                       &
    505509                THEN                                                           
    506510                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 )       
     
    512516                                 BTEST(wall_flags_static_0(k+1,j,i),1)     .AND.   &
    513517                                 BTEST(wall_flags_static_0(k_pp,j,i),1)    .AND.   &
    514                            .NOT. flag_set                           .OR.           &
    515                          k == nzt - 1 )                                            &
     518                           .NOT. flag_set                                  .OR.    &
     519                         ( k == nzt - 1 .AND. symmetry_flag == 0 ) )               &
    516520                THEN                                                           
    517521                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 )       
     
    635639                             BTEST(wall_flags_static_0(k+1,j,i),2)          .AND.     &
    636640                             BTEST(wall_flags_static_0(k,j,i),2) )          .OR.      &
    637                      k == nzt )                                                       &
     641                     ( k == nzt .AND. symmetry_flag == 0 ) )                          &
    638642                THEN                                                           
    639643                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 )     
     
    646650                                 BTEST(wall_flags_static_0(k_pp,j,i),2)    .AND.      &
    647651                           .NOT. flag_set                                  .OR.       &
    648                          k == nzt - 1 )                                               &
     652                         ( k == nzt - 1 .AND. symmetry_flag == 0 ) )                  &
    649653                THEN                                                           
    650654                   advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 )     
     
    10821086                                                                               
    10831087                flag_set = .FALSE.                                             
    1084                 IF ( .NOT. BTEST(wall_flags_static_0(k-1,j,i),0)  .AND.               &
    1085                            BTEST(wall_flags_static_0(k,j,i),0)    .OR.                &
    1086                      .NOT. BTEST(wall_flags_static_0(k_pp,j,i),0) .AND.               &                             
    1087                            BTEST(wall_flags_static_0(k,j,i),0)    .OR.                &
    1088                      k == nzt )                                                &
     1088                IF ( ( .NOT. BTEST(wall_flags_static_0(k-1,j,i),0)          .AND.     &
     1089                             BTEST(wall_flags_static_0(k,j,i),0)            .AND.     &
     1090                             BTEST(wall_flags_static_0(k+1,j,i),0) )        .OR.      &
     1091                     ( .NOT. BTEST(wall_flags_static_0(k_pp,j,i),0)         .AND.     &                             
     1092                             BTEST(wall_flags_static_0(k+1,j,i),0)          .AND.     &
     1093                             BTEST(wall_flags_static_0(k,j,i),0) )          .OR.      &
     1094                     ( k == nzt .AND. symmetry_flag == 0 ) )                          &
    10891095                THEN                                                           
    10901096                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 )             
     
    10961102                                 BTEST(wall_flags_static_0(k+1,j,i),0)     .AND.      &
    10971103                                 BTEST(wall_flags_static_0(k_pp,j,i),0)    .AND.      &
    1098                            .NOT. flag_set                           .OR.       &
    1099                          k == nzt - 1 )                                        &
     1104                           .NOT. flag_set                                  .OR.       &
     1105                         ( k == nzt - 1 .AND. symmetry_flag == 0 ) )                  &
    11001106                THEN                                                           
    11011107                   advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 )             
     
    15951601       ENDDO       
    15961602       
    1597        DO  k = nzt-1, nzt
     1603       DO  k = nzt-1, nzt-symmetry_flag
    15981604          ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
    15991605          ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp )
     
    16381644       
    16391645!
    1640 !--    Set resolved/turbulent flux at model top to zero (w-level)
     1646!--    Set resolved/turbulent flux at model top to zero (w-level). In case that
     1647!--    a symmetric behavior between bottom and top shall be guaranteed (closed
     1648!--    channel flow), the flux at nzt is also set to zero.
     1649       IF ( symmetry_flag == 1 ) THEN
     1650          flux_t(nzt) = 0.0_wp
     1651          diss_t(nzt) = 0.0_wp
     1652       ENDIF
    16411653       flux_t(nzt+1) = 0.0_wp
    16421654       diss_t(nzt+1) = 0.0_wp
     1655       
    16431656       
    16441657       IF ( limiter )  THEN
     
    23232336       ENDDO
    23242337       
    2325        DO  k = nzt-1, nzt
     2338       DO  k = nzt-1, nzt-symmetry_flag
    23262339!
    23272340!--       k index has to be modified near bottom and top, else array
     
    23682381       
    23692382!
    2370 !--    Set resolved/turbulent flux at model top to zero (w-level)
     2383!--    Set resolved/turbulent flux at model top to zero (w-level). In case that
     2384!--    a symmetric behavior between bottom and top shall be guaranteed (closed
     2385!--    channel flow), the flux at nzt is also set to zero.
     2386       IF ( symmetry_flag == 1 ) THEN
     2387          flux_t(nzt) = 0.0_wp
     2388          diss_t(nzt) = 0.0_wp
     2389          w_comp(nzt) = 0.0_wp
     2390       ENDIF
    23712391       flux_t(nzt+1) = 0.0_wp
    23722392       diss_t(nzt+1) = 0.0_wp
     
    28872907       ENDDO
    28882908       
    2889        DO  k = nzt-1, nzt
     2909       DO  k = nzt-1, nzt-symmetry_flag
    28902910!
    28912911!--       k index has to be modified near bottom and top, else array
     
    29322952       
    29332953!
    2934 !--    Set resolved/turbulent flux at model top to zero (w-level)
     2954!--    Set resolved/turbulent flux at model top to zero (w-level). In case that
     2955!--    a symmetric behavior between bottom and top shall be guaranteed (closed
     2956!--    channel flow), the flux at nzt is also set to zero.
     2957       IF ( symmetry_flag == 1 ) THEN
     2958          flux_t(nzt) = 0.0_wp
     2959          diss_t(nzt) = 0.0_wp
     2960          w_comp(nzt) = 0.0_wp
     2961       ENDIF
    29352962       flux_t(nzt+1) = 0.0_wp
    29362963       diss_t(nzt+1) = 0.0_wp
  • palm/trunk/SOURCE/check_parameters.f90

    r4331 r4340  
    2525! -----------------
    2626! 4172 2019-08-20 11:55:33Z oliver.maas
     27! Checks for closed channel flow implemented
     28!
     29! 11:55:33Z oliver.maas
    2730! Move 2-m potential temperature output to diagnostic_output_quantities
    2831!
     
    498501          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
    499502       ENDIF
    500        IF ( .NOT. constant_flux_layer )  THEN
     503       IF ( .NOT. constant_flux_layer .AND. topography /= 'closed_channel' )   &
     504       THEN
    501505          WRITE( action, '(A)' )  'constant_flux_layer = .FALSE.'
    502506       ENDIF
    503507       IF ( action /= ' ' )  THEN
    504           message_string = 'a non-flat topography does not allow ' //          &
     508          message_string = 'The specified topography does not allow ' //       &
    505509                           TRIM( action )
    506510          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
    507511       ENDIF
    508 
     512!
     513!--    Check illegal/untested parameter combinations for closed channel
     514       If ( topography == 'closed_channel' ) THEN
     515          symmetry_flag = 1
     516          message_string = 'Bottom and top boundary are treated equal'
     517          CALL message( 'check_parameters', 'PA0676', 0, 0, 0, 6, 0 )
     518       
     519          IF ( dz(1) /= dz(COUNT( dz /= -1.0_wp )) .OR.                        &
     520               dz_stretch_level /= -9999999.9_wp) THEN
     521             WRITE( message_string, * )  'dz should be equal close to the ' // &
     522                                         'boundaries due to symmetrical problem'
     523             CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
     524          ENDIF
     525       
     526          IF ( constant_flux_layer ) THEN
     527             WRITE( message_string, * )  'A constant flux layer is not '//     &
     528                                         'allowed if a closed channel '//      &
     529                                         'shall be used'
     530             CALL message( 'check_parameters', 'PA0688', 1, 2, 0, 6, 0 )
     531          ENDIF
     532       
     533          IF ( ocean_mode ) THEN
     534             WRITE( message_string, * )  'The ocean mode is not allowed if '// &
     535                                         'a closed channel shall be used'
     536             CALL message( 'check_parameters', 'PA0689', 1, 2, 0, 6, 0 )
     537          ENDIF
     538       
     539          IF ( momentum_advec /= 'ws-scheme' .OR.                              &
     540               scalar_advec /= 'ws-scheme' ) THEN
     541             WRITE( message_string, * )  'A closed channel require the '//     &
     542                                         'upwind scheme of Wicker and ' //     &
     543                                         'Skamarock as the advection scheme'
     544             CALL message( 'check_parameters', 'PA0690', 1, 2, 0, 6, 0 )
     545          ENDIF
     546       ENDIF
    509547    ENDIF
    510548
  • palm/trunk/SOURCE/init_grid.f90

    r4329 r4340  
    2525! -----------------
    2626! $Id$
     27! Topography closed channel flow with symmetric boundaries implemented
     28!
     29! 4329 2019-12-10 15:46:36Z motisi
    2730! Renamed wall_flags_0 to wall_flags_static_0
    2831!
     
    132135               momentum_advec, number_stretch_level_end,                       &
    133136               number_stretch_level_start, ocean_mode, psolver, scalar_advec,  &
    134                topography, use_surface_fluxes
     137               symmetry_flag, topography, use_surface_fluxes
    135138         
    136139    USE grid_variables,                                                        &
     
    360363!--    the first u/v- and w-level (k=0) are defined at same height (z=0).
    361364!--    The second u-level (k=1) corresponds to the top of the
    362 !--    Prandtl-layer.
    363        IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
     365!--    Prandtl-layer. In case of symmetric boundaries (closed channel flow),
     366!--    the first grid point is always at z=0.
     367       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 .OR.                              &
     368            topography == 'closed_channel' ) THEN
    364369          zu(0) = 0.0_wp
    365370       ELSE
     
    379384!--    The default value of dz_stretch_level_start is negative, thus the first
    380385!--    condition is always true. Hence, the second condition is necessary.
    381        DO  k = 2, nzt+1
     386       DO  k = 2, nzt+1-symmetry_flag
    382387          IF ( dz_stretch_level_start(n) <= zu(k-1) .AND.                      &
    383388               dz_stretch_level_start(n) /= -9999999.9_wp ) THEN
     
    408413          ENDIF
    409414       ENDDO
     415       
     416!
     417!--    If a closed channel flow is simulated, make sure that grid structure is 
     418!--    the same for both bottom and top boundary. (Hint: Using a different dz
     419!--    at the bottom and at the top makes no sense due to symmetric boundaries
     420!--    where dz should be equal. Therefore, different dz at the bottom and top 
     421!--    causes an abort (see check_parameters).)
     422       IF ( topography == 'closed_channel' ) THEN
     423          zu(nzt+1) = zu(nzt) + dz(1) * 0.5_wp
     424       ENDIF
    410425
    411426!
     
    413428!--    corresponding u-levels. In case of dirichlet bc for u and v at the
    414429!--    ground the first u- and w-level (k=0) are defined at same height (z=0).
    415 !--    The top w-level is extrapolated linearly.
     430!--    Per default, the top w-level is extrapolated linearly. In case of
     431!--    a closed channel flow, zu(nzt+1) and zw(nzt) must be set explicitely.
     432!--    (Hint: Using a different dz at the bottom and at the top makes no sense
     433!--    due to symmetric boundaries where dz should be equal. Therefore,
     434!--    different dz at the bottom and top causes an abort (see
     435!--    check_parameters).)
    416436       zw(0) = 0.0_wp
    417        DO  k = 1, nzt
     437       DO  k = 1, nzt-symmetry_flag
    418438          zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
    419439       ENDDO
    420        zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1) - zw(nzt) )
     440       IF ( topography == 'closed_channel' ) THEN
     441          zw(nzt)   = zw(nzt-1) + dz(1)
     442          zw(nzt+1) = zw(nzt) + dz(1)
     443       ELSE
     444          zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1) - zw(nzt) )
     445       ENDIF
    421446
    422447    ELSE
     
    18171842!   
    18181843!--       Initialilize 3D topography array, used later for initializing flags
    1819           topo(nzb+1:nzt+1,:,:) = IBSET( topo(nzb+1:nzt+1,:,:), 0 )
     1844          topo(nzb+1:nzt+1,:,:) = IBSET( topo(nzb+1:nzt+1,:,:), 0 )
     1845         
     1846       CASE ( 'closed_channel' )
     1847!   
     1848!--       Initialilize 3D topography array, used later for initializing flags
     1849          topo(nzb+1:nzt,:,:) = IBSET( topo(nzb+1:nzt,:,:), 0 )
    18201850
    18211851       CASE ( 'single_building' )
     
    22382268!--    is applicable. If this is not possible, abort.
    22392269       IF ( TRIM( topography_grid_convention ) == ' ' )  THEN
    2240           IF ( TRIM( topography ) /= 'single_building' .AND.                   &
     2270          IF ( TRIM( topography ) /= 'closed_channel' .AND.                    &
     2271               TRIM( topography ) /= 'single_building' .AND.                   &
    22412272               TRIM( topography ) /= 'single_street_canyon' .AND.              &
    22422273               TRIM( topography ) /= 'tunnel'  .AND.                           &
     
    22502281               'is not set. Its default value is & only valid for ',           &
    22512282               '"topography" = ''single_building'', ''tunnel'' ',              &
    2252                '''single_street_canyon'' & or ''read_from_file''.',            &
     2283               '''single_street_canyon'', ''closed_channel'' & or ',           &
     2284               '''read_from_file''.',                                          &
    22532285               '& Choose ''cell_edge'' or ''cell_center''.'
    22542286             CALL message( 'init_grid', 'PA0239', 1, 2, 0, 6, 0 )
     
    23442376!
    23452377!--          scalar grid
    2346              IF ( BTEST( topo(k,j,i), 0 ) )                                 &
     2378             IF ( BTEST( topo(k,j,i), 0 ) )                                    &
    23472379                wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 0 )
    23482380!
    23492381!--          u grid
    2350              IF ( BTEST( topo(k,j,i),   0 )  .AND.                          &
    2351                   BTEST( topo(k,j,i-1), 0 ) )                               &
     2382             IF ( BTEST( topo(k,j,i),   0 )  .AND.                             &
     2383                  BTEST( topo(k,j,i-1), 0 ) )                                  &
    23522384                wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 1 )
    23532385!
    23542386!--          v grid
    2355              IF ( BTEST( topo(k,j,i),   0 )  .AND.                          &
    2356                   BTEST( topo(k,j-1,i), 0 ) )                               &
     2387             IF ( BTEST( topo(k,j,i),   0 )  .AND.                             &
     2388                  BTEST( topo(k,j-1,i), 0 ) )                                  &
    23572389                 wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 2 )
    23582390
     
    23622394!
    23632395!--          w grid
    2364              IF ( BTEST( topo(k,j,i),   0 )  .AND.                          &
    2365                   BTEST( topo(k+1,j,i), 0 ) )                               &
     2396             IF ( BTEST( topo(k,j,i),   0 )  .AND.                             &
     2397                  BTEST( topo(k+1,j,i), 0 ) )                                  &
    23662398                wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 3 )
    23672399          ENDDO
    2368           wall_flags_static_0(nzt+1,j,i) = IBSET( wall_flags_static_0(nzt+1,j,i), 3 )
     2400         
     2401          IF ( topography /= 'closed_channel' ) THEN
     2402             wall_flags_static_0(nzt+1,j,i) = IBSET( wall_flags_static_0(nzt+1,j,i), 3 )
     2403          ENDIF
    23692404
    23702405       ENDDO
     
    24792514!
    24802515!--       Flags indicating downward facing walls
    2481           DO k = nzb+1, nzt
     2516          DO k = nzb+1, nzt+1
    24822517!
    24832518!--          Scalar grid
  • palm/trunk/SOURCE/modules.f90

    r4331 r4340  
    2525! -----------------
    2626! $Id$
     27! Flag for topography closed channel flow with symmetric boundaries introduced
     28!
     29! 4331 2019-12-10 18:25:02Z suehring
    2730! - do_output_at_2m, pt_2m_av
    2831!
     
    621624    INTEGER(iwp) ::  runnr = 0                         !< number of run in job chain
    622625    INTEGER(iwp) ::  subdomain_size                    !< number of grid points in (3d) subdomain including ghost points
     626    INTEGER(iwp) ::  symmetry_flag = 0                 !< flag for sterring the symmetric behavior of the bottom and top boundary
    623627    INTEGER(iwp) ::  terminate_coupled = 0             !< switch for steering termination in case of coupled runs
    624628    INTEGER(iwp) ::  terminate_coupled_remote = 0      !< switch for steering termination in case of coupled runs (condition of the remote model)
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r4337 r4340  
    2727! Former revisions:
    2828! -----------------
    29 !
     29! $Id$
    3030! Albedo indices for building_surface_pars are now declared as parameters to
    3131! prevent an error if the gfortran compiler with -Werror=unused-value is used
    32 !
    33 !
    34 ! $Id$
     32
     33! 4291 2019-11-11 12:36:54Z moh.hefny
    3534! Enabled RTM in case of biometeorology even if there is no vertical
    3635! surfaces or 3D vegetation in the domain
Note: See TracChangeset for help on using the changeset viewer.