Home


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

Changeset 767 for palm


Ignore:
Timestamp:
Oct 14, 2011 6:39:12 AM (13 years ago)
Author:
raasch
Message:

New:
---

Flow field initialization with given (e.g. measured) profiles. Profile data
for u-,v-velocity components + respective heights are given with new
inipar-parameters u_profile, v_profile, and uv_heights. Final profiles are
calculated from these given profiles by linear interpolation.
(check_parameters, header, init_3d_model, modules, parin)

Changed:

ug,vg replaced by u_init,v_init as the Dirichlet top boundary condition
(boundary_conds)

dirichlet_0 conditions moved from init_3d_model to
check_parameters (check_parameters, init_3d_model)

Errors:

bugfix: dirichlet_0 conditions moved from init_3d_model to
check_parameters (check_parameters, init_3d_model)

Location:
palm/trunk/SOURCE
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/boundary_conds.f90

    r668 r767  
    44! Current revisions:
    55! -----------------
    6 !
     6! ug,vg replaced by u_init,v_init as the Dirichlet top boundary condition
    77! 
    88! Former revisions:
    … …  
    9090!--    Top boundary
    9191       IF ( ibc_uv_t == 0 )  THEN
    92            u_p(nzt+1,:,:) = ug(nzt+1)
    93            v_p(nzt+1,:,:) = vg(nzt+1)
     92           u_p(nzt+1,:,:) = u_init(nzt+1)
     93           v_p(nzt+1,:,:) = v_init(nzt+1)
    9494       ELSE
    9595           u_p(nzt+1,:,:) = u_p(nzt,:,:)
    … …  
    312312!--    Top boundary at the outflow
    313313       IF ( ibc_uv_t == 0 )  THEN
    314           u_p(nzt+1,-1,:) = ug(nzt+1)
    315           v_p(nzt+1,0,:)  = vg(nzt+1)
     314          u_p(nzt+1,-1,:) = u_init(nzt+1)
     315          v_p(nzt+1,0,:)  = v_init(nzt+1)
    316316       ELSE
    317317          u_p(nzt+1,-1,:) = u(nzt,-1,:)
    … …  
    410410!--    Top boundary at the outflow
    411411       IF ( ibc_uv_t == 0 )  THEN
    412           u_p(nzt+1,ny+1,:) = ug(nzt+1)
    413           v_p(nzt+1,ny+1,:) = vg(nzt+1)
     412          u_p(nzt+1,ny+1,:) = u_init(nzt+1)
     413          v_p(nzt+1,ny+1,:) = v_init(nzt+1)
    414414       ELSE
    415415          u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
    … …  
    508508!--    Top boundary at the outflow
    509509       IF ( ibc_uv_t == 0 )  THEN
    510           u_p(nzt+1,:,-1) = ug(nzt+1)
    511           v_p(nzt+1,:,-1) = vg(nzt+1)
     510          u_p(nzt+1,:,-1) = u_init(nzt+1)
     511          v_p(nzt+1,:,-1) = v_init(nzt+1)
    512512       ELSE
    513513          u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
    … …  
    606606!--    Top boundary at the outflow
    607607       IF ( ibc_uv_t == 0 )  THEN
    608           u_p(nzt+1,:,nx+1) = ug(nzt+1)
    609           v_p(nzt+1,:,nx+1) = vg(nzt+1)
     608          u_p(nzt+1,:,nx+1) = u_init(nzt+1)
     609          v_p(nzt+1,:,nx+1) = v_init(nzt+1)
    610610       ELSE
    611611          u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)

  • palm/trunk/SOURCE/check_parameters.f90

    r708 r767  
    44! Current revisions:
    55! -----------------
    6 !
     6! Calculating u,v-profiles from given profiles by linear interpolation.
     7! bugfix: dirichlet_0 conditions for ug/vg moved from init_3d_model to here
    78!
    89! Former revisions:
    … …  
    176177    CHARACTER (LEN=100) ::  action
    177178
    178     INTEGER ::  i, ilen, intervals, iremote = 0, iter, j, k, nnxh, nnyh, &
    179          position, prec
     179    INTEGER ::  i, ilen, intervals, iremote = 0, iter, j, k, kk, nnxh, nnyh, &
     180                position, prec
    180181    LOGICAL ::  found, ldum
    181182    REAL    ::  gradient, maxn, maxp, remote = 0.0, &
    … …  
    771772
    772773!
    773 !--    Initial profiles for 1D and 3D model, respectively
    774        u_init  = ug_surface
    775        v_init  = vg_surface
     774!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
    776775       pt_init = pt_surface
    777776       IF ( humidity )        q_init  = q_surface
    … …  
    838837       ENDIF
    839838
    840        u_init = ug
    841 
    842 !
    843 !--    In case of no given gradients for ug, choose a vanishing gradient
     839!
     840!--    In case of no given gradients for ug, choose a zero gradient
    844841       IF ( ug_vertical_gradient_level(1) == -9999999.9 )  THEN
    845842          ug_vertical_gradient_level(1) = 0.0
    … …  
    904901       ENDIF
    905902
    906        v_init = vg
    907  
    908 !
    909 !--    In case of no given gradients for vg, choose a vanishing gradient
     903!
     904!--    In case of no given gradients for vg, choose a zero gradient
    910905       IF ( vg_vertical_gradient_level(1) == -9999999.9 )  THEN
    911906          vg_vertical_gradient_level(1) = 0.0
     907       ENDIF
     908
     909!
     910!--    Let the initial wind profiles be the calculated ug/vg profiles or
     911!--    interpolate them from wind profile data (if given)
     912       IF ( u_profile(1) == 9999999.9  .AND.  v_profile(1) == 9999999.9 )  THEN
     913
     914          u_init = ug
     915          v_init = vg
     916
     917       ELSEIF ( u_profile(1) == 0.0  .AND.  v_profile(1) == 0.0 )  THEN
     918
     919          IF ( uv_heights(1) /= 0.0 )  THEN
     920             message_string = 'uv_heights(1) must be 0.0'
     921             CALL message( 'check_parameters', 'PA0345', 1, 2, 0, 6, 0 )
     922          ENDIF
     923
     924          use_prescribed_profile_data = .TRUE.
     925
     926          kk = 1
     927          u_init(0) = 0.0
     928          v_init(0) = 0.0
     929
     930          DO  k = 1, nz+1
     931
     932             IF ( kk < 100 )  THEN
     933                DO WHILE ( uv_heights(kk+1) <= zu(k) )
     934                   kk = kk + 1
     935                   IF ( kk == 100 )  EXIT
     936                ENDDO
     937             ENDIF
     938
     939             IF ( kk < 100 )  THEN
     940                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
     941                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
     942                                       ( u_profile(kk+1) - u_profile(kk) )
     943                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
     944                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
     945                                       ( v_profile(kk+1) - v_profile(kk) )
     946             ELSE
     947                u_init(k) = u_profile(kk)
     948                v_init(k) = v_profile(kk)
     949             ENDIF
     950
     951          ENDDO
     952
     953       ELSE
     954
     955          message_string = 'u_profile(1) and v_profile(1) must be 0.0'
     956          CALL message( 'check_parameters', 'PA0346', 1, 2, 0, 6, 0 )
     957
    912958       ENDIF
    913959
    … …  
    15391585       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
    15401586          ibc_uv_t = 0
     1587          IF ( bc_uv_t == 'dirichlet_0' )  THEN
     1588!
     1589!--          Velocities for the initial u,v-profiles are set zero at the top
     1590!--          in case of dirichlet_0 conditions
     1591             u_init(nzt+1)    = 0.0
     1592             v_init(nzt+1)    = 0.0
     1593          ENDIF
    15411594       ELSEIF ( bc_uv_t == 'neumann' )  THEN
    15421595          ibc_uv_t = 1

  • palm/trunk/SOURCE/header.f90

    r760 r767  
    44! Current revisions:
    55! -----------------
    6 !
     6! Output of given initial u,v-profiles
    77!
    88! Former revisions:
    … …  
    13171317    WRITE ( io, 424 )  TRIM( coordinates ), TRIM( vgcomponent ), &
    13181318                       TRIM( gradients ), TRIM( slices )
     1319
     1320!
     1321!-- Initial wind profiles
     1322    IF ( u_profile(1) /= 9999999.9 )  WRITE ( io, 427 )
    13191323
    13201324!
    … …  
    19251929            '       Gradient:    ',A,'  (m/s)/100m'/ &
    19261930            '       Gridpoint:   ',A)
     1931427 FORMAT (/'    Initial wind profiles (u,v) are interpolated from given'// &
     1932                  ' profiles')
    19271933450 FORMAT (//' LES / Turbulence quantities:'/ &
    19281934              ' ---------------------------'/)

  • palm/trunk/SOURCE/init_3d_model.f90

    r760 r767  
    77! Current revisions:
    88! ------------------
    9 !
     9! adjustments for prescribed u,v-profiles
     10! bugfix: dirichlet_0 conditions for ug/vg moved to check_parameters
    1011!
    1112! Former revisions:
    … …  
    694695!--    Apply channel flow boundary condition
    695696       IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
    696 
    697697          u(nzt+1,:,:) = 0.0
    698698          v(nzt+1,:,:) = 0.0
    699 
    700 !--       For the Dirichlet condition to be correctly applied at the top, set
    701 !--       ug and vg to zero there
    702           ug(nzt+1)    = 0.0
    703           vg(nzt+1)    = 0.0
    704 
    705699       ENDIF
    706700
    … …  
    937931    THEN
    938932!
    939 !--    When reading data for cyclic fill of 3D prerun data, first read
    940 !--    some of the global variables from restart file
     933!--    When reading data for cyclic fill of 3D prerun data files, read
     934!--    some of the global variables from the restart file which are required
     935!--    for initializing the inflow
    941936       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
    942937
    … …  
    951946          ENDDO
    952947
    953 !
    954 !--       Initialization of the turbulence recycling method
    955           IF ( turbulent_inflow )  THEN
    956 !
    957 !--          Store temporally and horizontally averaged vertical profiles to be
    958 !--          used as mean inflow profiles
    959              ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
    960 
    961              mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
    962              mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
    963              mean_inflow_profiles(:,4) = hom_sum(:,4,0)    ! pt
    964              mean_inflow_profiles(:,5) = hom_sum(:,8,0)    ! e
    965 
    966 !
    967 !--          Use these mean profiles for the inflow (provided that Dirichlet
    968 !--          conditions are used)
    969              IF ( inflow_l )  THEN
    970                 DO  j = nysg, nyng
    971                    DO  k = nzb, nzt+1
    972                       u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
    973                       v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
    974                       w(k,j,nxlg:-1)  = 0.0
    975                       pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
    976                       e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
    977                    ENDDO
    978                 ENDDO
    979              ENDIF
    980 
    981 !
    982 !--          Calculate the damping factors to be used at the inflow. For a
    983 !--          turbulent inflow the turbulent fluctuations have to be limited
    984 !--          vertically because otherwise the turbulent inflow layer will grow
    985 !--          in time.
    986              IF ( inflow_damping_height == 9999999.9 )  THEN
    987 !
    988 !--             Default: use the inversion height calculated by the prerun; if
    989 !--             this is zero, inflow_damping_height must be explicitly
    990 !--             specified.
    991                 IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0 )  THEN
    992                    inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
    993                 ELSE
    994                    WRITE( message_string, * ) 'inflow_damping_height must be ',&
    995                         'explicitly specified because&the inversion height ', &
    996                         'calculated by the prerun is zero.'
    997                    CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
    998                 ENDIF
    999 
    1000              ENDIF
    1001 
    1002              IF ( inflow_damping_width == 9999999.9 )  THEN
    1003 !
    1004 !--             Default for the transition range: one tenth of the undamped
    1005 !--             layer
    1006                 inflow_damping_width = 0.1 * inflow_damping_height
    1007 
    1008              ENDIF
    1009 
    1010              ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
    1011 
    1012              DO  k = nzb, nzt+1
    1013 
    1014                 IF ( zu(k) <= inflow_damping_height )  THEN
    1015                    inflow_damping_factor(k) = 1.0
    1016                 ELSEIF ( zu(k) <= inflow_damping_height +  &
    1017                                   inflow_damping_width )  THEN
    1018                    inflow_damping_factor(k) = 1.0 -                            &
    1019                                            ( zu(k) - inflow_damping_height ) / &
    1020                                            inflow_damping_width
    1021                 ELSE
    1022                    inflow_damping_factor(k) = 0.0
    1023                 ENDIF
    1024 
    1025              ENDDO
    1026           ENDIF
    1027 
    1028948       ENDIF
    1029949
    … …  
    1038958#endif
    1039959       ENDDO
     960
     961!
     962!--    Initialization of the turbulence recycling method
     963       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.  &
     964            turbulent_inflow )  THEN
     965!
     966!--       First store the profiles to be used at the inflow.
     967!--       These profiles are the (temporally) and horizontally averaged vertical
     968!--       profiles from the prerun. Alternatively, prescribed profiles
     969!--       for u,v-components can be used.
     970          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
     971
     972          IF ( use_prescribed_profile_data )  THEN
     973             mean_inflow_profiles(:,1) = u_init            ! u
     974             mean_inflow_profiles(:,2) = v_init            ! v
     975          ELSE
     976             mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
     977             mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
     978          ENDIF
     979          mean_inflow_profiles(:,4) = hom_sum(:,4,0)       ! pt
     980          mean_inflow_profiles(:,5) = hom_sum(:,8,0)       ! e
     981
     982!
     983!--       If necessary, adjust the horizontal flow field to the prescribed
     984!--       profiles
     985          IF ( use_prescribed_profile_data )  THEN
     986             DO  i = nxlg, nxrg
     987                DO  j = nysg, nyng
     988                   DO  k = nzb, nzt+1
     989                      u(k,j,i) = u(k,j,i) - hom_sum(k,1,0) + u_init(k)
     990                      v(k,j,i) = v(k,j,i) - hom_sum(k,2,0) + v_init(k)
     991                   ENDDO
     992                ENDDO
     993             ENDDO
     994          ENDIF
     995
     996!
     997!--       Use these mean profiles at the inflow (provided that Dirichlet
     998!--       conditions are used)
     999          IF ( inflow_l )  THEN
     1000             DO  j = nysg, nyng
     1001                DO  k = nzb, nzt+1
     1002                   u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
     1003                   v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
     1004                   w(k,j,nxlg:-1)  = 0.0
     1005                   pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
     1006                   e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
     1007                ENDDO
     1008             ENDDO
     1009          ENDIF
     1010
     1011!
     1012!--       Calculate the damping factors to be used at the inflow. For a
     1013!--       turbulent inflow the turbulent fluctuations have to be limited
     1014!--       vertically because otherwise the turbulent inflow layer will grow
     1015!--       in time.
     1016          IF ( inflow_damping_height == 9999999.9 )  THEN
     1017!
     1018!--          Default: use the inversion height calculated by the prerun; if
     1019!--          this is zero, inflow_damping_height must be explicitly
     1020!--          specified.
     1021             IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0 )  THEN
     1022                inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
     1023             ELSE
     1024                WRITE( message_string, * ) 'inflow_damping_height must be ',&
     1025                     'explicitly specified because&the inversion height ', &
     1026                     'calculated by the prerun is zero.'
     1027                CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
     1028             ENDIF
     1029
     1030          ENDIF
     1031
     1032          IF ( inflow_damping_width == 9999999.9 )  THEN
     1033!
     1034!--          Default for the transition range: one tenth of the undamped
     1035!--          layer
     1036             inflow_damping_width = 0.1 * inflow_damping_height
     1037
     1038          ENDIF
     1039
     1040          ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
     1041
     1042          DO  k = nzb, nzt+1
     1043
     1044             IF ( zu(k) <= inflow_damping_height )  THEN
     1045                inflow_damping_factor(k) = 1.0
     1046             ELSEIF ( zu(k) <= inflow_damping_height +  &
     1047                               inflow_damping_width )  THEN
     1048                inflow_damping_factor(k) = 1.0 -                            &
     1049                                        ( zu(k) - inflow_damping_height ) / &
     1050                                        inflow_damping_width
     1051             ELSE
     1052                inflow_damping_factor(k) = 0.0
     1053             ENDIF
     1054
     1055          ENDDO
     1056
     1057       ENDIF
    10401058
    10411059!
    … …  
    11461164    IF ( conserve_volume_flow )  THEN
    11471165
    1148        IF  ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
     1166       IF ( use_prescribed_profile_data )  THEN
     1167
     1168          volume_flow_initial_l = 0.0
     1169          volume_flow_area_l    = 0.0
     1170
     1171          IF ( nxr == nx )  THEN
     1172             DO  j = nys, nyn
     1173                DO  k = nzb_2d(j,nx)+1, nzt
     1174                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
     1175                                              u_init(k) * dzw(k)
     1176                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
     1177                ENDDO
     1178             ENDDO
     1179          ENDIF
     1180         
     1181          IF ( nyn == ny )  THEN
     1182             DO  i = nxl, nxr
     1183                DO  k = nzb_2d(ny,i)+1, nzt 
     1184                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
     1185                                              v_init(k) * dzw(k)
     1186                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
     1187                ENDDO
     1188             ENDDO
     1189          ENDIF
     1190
     1191#if defined( __parallel )
     1192          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
     1193                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
     1194          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
     1195                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
     1196
     1197#else
     1198          volume_flow_initial = volume_flow_initial_l
     1199          volume_flow_area    = volume_flow_area_l
     1200#endif 
     1201
     1202       ELSEIF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
    11491203
    11501204          volume_flow_initial_l = 0.0

  • palm/trunk/SOURCE/modules.f90

    r760 r767  
    55! Current revisions:
    66! -----------------
    7 !
     7! +u_profile, v_profile, uv_heights, use_prescribed_profile_data
    88!
    99! Former revisions:
    … …  
    536536                stop_dt = .FALSE., synchronous_exchange = .FALSE., &
    537537                terminate_run = .FALSE., turbulent_inflow = .FALSE., &
     538                use_prescribed_profile_data = .FALSE., &
    538539                use_prior_plot1d_parameters = .FALSE., use_reference = .FALSE.,&
    539540                use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
    540541                use_ug_for_galilei_tr = .TRUE., use_upstream_for_tke = .FALSE.,&
    541                 wall_adjustment = .TRUE., ws_scheme_sca = .FALSE.,             &
     542                wall_adjustment = .TRUE., ws_scheme_sca = .FALSE., &
    542543                ws_scheme_mom = .FALSE.
    543544
    … …  
    636637             time_domask(max_masks) = 0.0, &
    637638             tsc(10) = (/ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /), &
     639             u_profile(100) = 9999999.9, uv_heights(100) = 9999999.9, &
     640             v_profile(100) = 9999999.9, &
    638641             ug_vertical_gradient(10) = 0.0, &
    639642             ug_vertical_gradient_level(10) = -9999999.9, &

  • palm/trunk/SOURCE/parin.f90

    r760 r767  
    44! Current revisions:
    55! -----------------
    6 !
     6! +u_profile, v_profile, uv_heights in inipar
    77!
    88! Former revisions:
    … …  
    166166             top_momentumflux_u, top_momentumflux_v, top_salinityflux, &
    167167             turbulent_inflow, ug_surface, ug_vertical_gradient, &
    168              ug_vertical_gradient_level, u_bulk, ups_limit_e, ups_limit_pt, &
     168             ug_vertical_gradient_level, ups_limit_e, ups_limit_pt, &
    169169             ups_limit_u, ups_limit_v, ups_limit_w, use_surface_fluxes, &
    170170             use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke, &
    171              vg_surface, vg_vertical_gradient, vg_vertical_gradient_level, &
    172              v_bulk, wall_adjustment, wall_heatflux, wall_humidityflux, &
    173              wall_scalarflux
     171             uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient, &
     172             vg_vertical_gradient_level, v_bulk, v_profile, wall_adjustment, &
     173             wall_heatflux, wall_humidityflux, wall_scalarflux
    174174
    175175
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |