Changeset 767 for palm/trunk/SOURCE
- Timestamp:
- Oct 14, 2011 6:39:12 AM (13 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 6 edited
-
boundary_conds.f90 (modified) (6 diffs)
-
check_parameters.f90 (modified) (6 diffs)
-
header.f90 (modified) (3 diffs)
-
init_3d_model.f90 (modified) (6 diffs)
-
modules.f90 (modified) (3 diffs)
-
parin.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/boundary_conds.f90
r668 r767 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! ug,vg replaced by u_init,v_init as the Dirichlet top boundary condition 7 7 ! 8 8 ! Former revisions: … … 90 90 !-- Top boundary 91 91 IF ( ibc_uv_t == 0 ) THEN 92 u_p(nzt+1,:,:) = u g(nzt+1)93 v_p(nzt+1,:,:) = v g(nzt+1)92 u_p(nzt+1,:,:) = u_init(nzt+1) 93 v_p(nzt+1,:,:) = v_init(nzt+1) 94 94 ELSE 95 95 u_p(nzt+1,:,:) = u_p(nzt,:,:) … … 312 312 !-- Top boundary at the outflow 313 313 IF ( ibc_uv_t == 0 ) THEN 314 u_p(nzt+1,-1,:) = u g(nzt+1)315 v_p(nzt+1,0,:) = v g(nzt+1)314 u_p(nzt+1,-1,:) = u_init(nzt+1) 315 v_p(nzt+1,0,:) = v_init(nzt+1) 316 316 ELSE 317 317 u_p(nzt+1,-1,:) = u(nzt,-1,:) … … 410 410 !-- Top boundary at the outflow 411 411 IF ( ibc_uv_t == 0 ) THEN 412 u_p(nzt+1,ny+1,:) = u g(nzt+1)413 v_p(nzt+1,ny+1,:) = v g(nzt+1)412 u_p(nzt+1,ny+1,:) = u_init(nzt+1) 413 v_p(nzt+1,ny+1,:) = v_init(nzt+1) 414 414 ELSE 415 415 u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:) … … 508 508 !-- Top boundary at the outflow 509 509 IF ( ibc_uv_t == 0 ) THEN 510 u_p(nzt+1,:,-1) = u g(nzt+1)511 v_p(nzt+1,:,-1) = v g(nzt+1)510 u_p(nzt+1,:,-1) = u_init(nzt+1) 511 v_p(nzt+1,:,-1) = v_init(nzt+1) 512 512 ELSE 513 513 u_p(nzt+1,:,-1) = u_p(nzt,:,-1) … … 606 606 !-- Top boundary at the outflow 607 607 IF ( ibc_uv_t == 0 ) THEN 608 u_p(nzt+1,:,nx+1) = u g(nzt+1)609 v_p(nzt+1,:,nx+1) = v g(nzt+1)608 u_p(nzt+1,:,nx+1) = u_init(nzt+1) 609 v_p(nzt+1,:,nx+1) = v_init(nzt+1) 610 610 ELSE 611 611 u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1) -
palm/trunk/SOURCE/check_parameters.f90
r708 r767 4 4 ! Current revisions: 5 5 ! ----------------- 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 7 8 ! 8 9 ! Former revisions: … … 176 177 CHARACTER (LEN=100) :: action 177 178 178 INTEGER :: i, ilen, intervals, iremote = 0, iter, j, k, nnxh, nnyh, &179 position, prec179 INTEGER :: i, ilen, intervals, iremote = 0, iter, j, k, kk, nnxh, nnyh, & 180 position, prec 180 181 LOGICAL :: found, ldum 181 182 REAL :: gradient, maxn, maxp, remote = 0.0, & … … 771 772 772 773 ! 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) 776 775 pt_init = pt_surface 777 776 IF ( humidity ) q_init = q_surface … … 838 837 ENDIF 839 838 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 844 841 IF ( ug_vertical_gradient_level(1) == -9999999.9 ) THEN 845 842 ug_vertical_gradient_level(1) = 0.0 … … 904 901 ENDIF 905 902 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 910 905 IF ( vg_vertical_gradient_level(1) == -9999999.9 ) THEN 911 906 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 912 958 ENDIF 913 959 … … 1539 1585 IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' ) THEN 1540 1586 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 1541 1594 ELSEIF ( bc_uv_t == 'neumann' ) THEN 1542 1595 ibc_uv_t = 1 -
palm/trunk/SOURCE/header.f90
r760 r767 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Output of given initial u,v-profiles 7 7 ! 8 8 ! Former revisions: … … 1317 1317 WRITE ( io, 424 ) TRIM( coordinates ), TRIM( vgcomponent ), & 1318 1318 TRIM( gradients ), TRIM( slices ) 1319 1320 ! 1321 !-- Initial wind profiles 1322 IF ( u_profile(1) /= 9999999.9 ) WRITE ( io, 427 ) 1319 1323 1320 1324 ! … … 1925 1929 ' Gradient: ',A,' (m/s)/100m'/ & 1926 1930 ' Gridpoint: ',A) 1931 427 FORMAT (/' Initial wind profiles (u,v) are interpolated from given'// & 1932 ' profiles') 1927 1933 450 FORMAT (//' LES / Turbulence quantities:'/ & 1928 1934 ' ---------------------------'/) -
palm/trunk/SOURCE/init_3d_model.f90
r760 r767 7 7 ! Current revisions: 8 8 ! ------------------ 9 ! 9 ! adjustments for prescribed u,v-profiles 10 ! bugfix: dirichlet_0 conditions for ug/vg moved to check_parameters 10 11 ! 11 12 ! Former revisions: … … 694 695 !-- Apply channel flow boundary condition 695 696 IF ( TRIM( bc_uv_t ) == 'dirichlet_0' ) THEN 696 697 697 u(nzt+1,:,:) = 0.0 698 698 v(nzt+1,:,:) = 0.0 699 700 !-- For the Dirichlet condition to be correctly applied at the top, set701 !-- ug and vg to zero there702 ug(nzt+1) = 0.0703 vg(nzt+1) = 0.0704 705 699 ENDIF 706 700 … … 937 931 THEN 938 932 ! 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 941 936 IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN 942 937 … … 951 946 ENDDO 952 947 953 !954 !-- Initialization of the turbulence recycling method955 IF ( turbulent_inflow ) THEN956 !957 !-- Store temporally and horizontally averaged vertical profiles to be958 !-- used as mean inflow profiles959 ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )960 961 mean_inflow_profiles(:,1) = hom_sum(:,1,0) ! u962 mean_inflow_profiles(:,2) = hom_sum(:,2,0) ! v963 mean_inflow_profiles(:,4) = hom_sum(:,4,0) ! pt964 mean_inflow_profiles(:,5) = hom_sum(:,8,0) ! e965 966 !967 !-- Use these mean profiles for the inflow (provided that Dirichlet968 !-- conditions are used)969 IF ( inflow_l ) THEN970 DO j = nysg, nyng971 DO k = nzb, nzt+1972 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.0975 pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)976 e(k,j,nxlg:-1) = mean_inflow_profiles(k,5)977 ENDDO978 ENDDO979 ENDIF980 981 !982 !-- Calculate the damping factors to be used at the inflow. For a983 !-- turbulent inflow the turbulent fluctuations have to be limited984 !-- vertically because otherwise the turbulent inflow layer will grow985 !-- in time.986 IF ( inflow_damping_height == 9999999.9 ) THEN987 !988 !-- Default: use the inversion height calculated by the prerun; if989 !-- this is zero, inflow_damping_height must be explicitly990 !-- specified.991 IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0 ) THEN992 inflow_damping_height = hom_sum(nzb+6,pr_palm,0)993 ELSE994 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 ENDIF999 1000 ENDIF1001 1002 IF ( inflow_damping_width == 9999999.9 ) THEN1003 !1004 !-- Default for the transition range: one tenth of the undamped1005 !-- layer1006 inflow_damping_width = 0.1 * inflow_damping_height1007 1008 ENDIF1009 1010 ALLOCATE( inflow_damping_factor(nzb:nzt+1) )1011 1012 DO k = nzb, nzt+11013 1014 IF ( zu(k) <= inflow_damping_height ) THEN1015 inflow_damping_factor(k) = 1.01016 ELSEIF ( zu(k) <= inflow_damping_height + &1017 inflow_damping_width ) THEN1018 inflow_damping_factor(k) = 1.0 - &1019 ( zu(k) - inflow_damping_height ) / &1020 inflow_damping_width1021 ELSE1022 inflow_damping_factor(k) = 0.01023 ENDIF1024 1025 ENDDO1026 ENDIF1027 1028 948 ENDIF 1029 949 … … 1038 958 #endif 1039 959 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 1040 1058 1041 1059 ! … … 1146 1164 IF ( conserve_volume_flow ) THEN 1147 1165 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 1149 1203 1150 1204 volume_flow_initial_l = 0.0 -
palm/trunk/SOURCE/modules.f90
r760 r767 5 5 ! Current revisions: 6 6 ! ----------------- 7 ! 7 ! +u_profile, v_profile, uv_heights, use_prescribed_profile_data 8 8 ! 9 9 ! Former revisions: … … 536 536 stop_dt = .FALSE., synchronous_exchange = .FALSE., & 537 537 terminate_run = .FALSE., turbulent_inflow = .FALSE., & 538 use_prescribed_profile_data = .FALSE., & 538 539 use_prior_plot1d_parameters = .FALSE., use_reference = .FALSE.,& 539 540 use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., & 540 541 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., & 542 543 ws_scheme_mom = .FALSE. 543 544 … … 636 637 time_domask(max_masks) = 0.0, & 637 638 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, & 638 641 ug_vertical_gradient(10) = 0.0, & 639 642 ug_vertical_gradient_level(10) = -9999999.9, & -
palm/trunk/SOURCE/parin.f90
r760 r767 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +u_profile, v_profile, uv_heights in inipar 7 7 ! 8 8 ! Former revisions: … … 166 166 top_momentumflux_u, top_momentumflux_v, top_salinityflux, & 167 167 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, & 169 169 ups_limit_u, ups_limit_v, ups_limit_w, use_surface_fluxes, & 170 170 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_ scalarflux171 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 174 174 175 175
Note: See TracChangeset
for help on using the changeset viewer.
