Changeset 4281 for palm/trunk
- Timestamp:
- Oct 29, 2019 3:15:39 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 1 deleted
- 4 edited
-
Makefile (modified) (4 diffs)
-
boundary_conds.f90 (deleted)
-
dynamics_mod.f90 (modified) (5 diffs)
-
module_interface.f90 (modified) (4 diffs)
-
time_integration.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r4270 r4281 25 25 # ----------------- 26 26 # $Id$ 27 # delete boundary_conds, added missing dependencies 28 # 29 # 4270 2019-10-23 10:46:20Z monakurppa 27 30 # - Implement offline nesting for salsa and add dependency of nesting_offl_mod 28 31 # on salsa_mod … … 127 130 basic_constants_and_equations_mod.f90 \ 128 131 biometeorology_mod.f90 \ 129 boundary_conds.f90 \130 132 buoyancy.f90 \ 131 133 calc_mean_profile.f90 \ … … 345 347 palm_date_time_mod.o \ 346 348 radiation_model_mod.o 347 boundary_conds.o: \348 bulk_cloud_model_mod.o \349 chemistry_model_mod.o \350 mod_kinds.o \351 modules.o \352 salsa_mod.o \353 surface_mod.o \354 turbulence_closure_mod.o355 349 bulk_cloud_model_mod.o: \ 356 350 basic_constants_and_equations_mod.o \ … … 553 547 dynamics_mod.o: \ 554 548 mod_kinds.o \ 549 surface_mod.o \ 550 pmc_interface_mod.o \ 555 551 modules.o 556 552 exchange_horiz.o: \ -
palm/trunk/SOURCE/dynamics_mod.f90
r4097 r4281 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Moved boundary conditions in dynamics module 28 ! 29 ! 4097 2019-07-15 11:59:11Z suehring 27 30 ! Avoid overlong lines - limit is 132 characters per line 28 31 ! … … 39 42 40 43 USE arrays_3d, & 41 ONLY: pt, pt_1, pt_2, pt_p, & 44 ONLY: c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l, & 45 dzu, & 46 pt, pt_1, pt_2, pt_init, pt_p, & 42 47 q, q_1, q_2, q_p, & 43 48 s, s_1, s_2, s_p, & 44 u, u_1, u_2, u_ p, &45 v, v_1, v_2, v_p, &46 w, w_1, w_2, w_p 49 u, u_1, u_2, u_init, u_p, u_m_l, u_m_n, u_m_r, u_m_s, & 50 v, v_1, v_2, v_p, v_init, v_m_l, v_m_n, v_m_r, v_m_s, & 51 w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s 47 52 48 53 USE control_parameters, & 49 ONLY: length, & 54 ONLY: bc_dirichlet_l, & 55 bc_dirichlet_s, & 56 bc_radiation_l, & 57 bc_radiation_n, & 58 bc_radiation_r, & 59 bc_radiation_s, & 60 bc_pt_t_val, & 61 bc_q_t_val, & 62 bc_s_t_val, & 63 child_domain, & 64 coupling_mode, & 65 dt_3d, & 66 ibc_pt_b, & 67 ibc_pt_t, & 68 ibc_q_b, & 69 ibc_q_t, & 70 ibc_s_b, & 71 ibc_s_t, & 72 ibc_uv_b, & 73 ibc_uv_t, & 74 intermediate_timestep_count, & 75 length, & 76 nesting_offline, & 77 nudging, & 50 78 restart_string, & 51 79 humidity, & 52 80 neutral, & 53 passive_scalar 81 passive_scalar, & 82 tsc, & 83 use_cmax 84 85 USE grid_variables, & 86 ONLY: ddx, & 87 ddy, & 88 dx, & 89 dy 54 90 55 91 USE indices, & 56 92 ONLY: nbgp, & 93 nx, & 57 94 nxl, & 95 nxlg, & 58 96 nxr, & 97 nxrg, & 98 ny, & 59 99 nys, & 100 nysg, & 60 101 nyn, & 102 nyng, & 61 103 nzb, & 62 104 nzt 63 105 64 106 USE kinds 107 108 USE pegrid 109 110 USE pmc_interface, & 111 ONLY : nesting_mode 112 113 USE surface_mod, & 114 ONLY : bc_h 115 65 116 66 117 IMPLICIT NONE … … 90 141 dynamics_exchange_horiz, & 91 142 dynamics_prognostic_equations, & 143 dynamics_boundary_conditions, & 92 144 dynamics_swap_timelevel, & 93 145 dynamics_3d_data_averaging, & … … 169 221 END INTERFACE dynamics_prognostic_equations 170 222 223 INTERFACE dynamics_boundary_conditions 224 MODULE PROCEDURE dynamics_boundary_conditions 225 END INTERFACE dynamics_boundary_conditions 226 171 227 INTERFACE dynamics_swap_timelevel 172 228 MODULE PROCEDURE dynamics_swap_timelevel … … 657 713 658 714 715 !--------------------------------------------------------------------------------------------------! 716 ! Description: 717 ! ------------ 718 !> Compute boundary conditions of dynamics model 719 !--------------------------------------------------------------------------------------------------! 720 SUBROUTINE dynamics_boundary_conditions 721 722 IMPLICIT NONE 723 724 INTEGER(iwp) :: i !< grid index x direction 725 INTEGER(iwp) :: j !< grid index y direction 726 INTEGER(iwp) :: k !< grid index z direction 727 INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls 728 INTEGER(iwp) :: m !< running index surface elements 729 730 REAL(wp) :: c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition 731 REAL(wp) :: denom !< horizontal gradient of velocity component normal to the outflow boundary 732 733 ! 734 !-- Bottom boundary 735 IF ( ibc_uv_b == 1 ) THEN 736 u_p(nzb,:,:) = u_p(nzb+1,:,:) 737 v_p(nzb,:,:) = v_p(nzb+1,:,:) 738 ENDIF 739 ! 740 !-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case 741 !-- of downward-facing surfaces. 742 DO l = 0, 1 743 !$OMP PARALLEL DO PRIVATE( i, j, k ) 744 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 745 !$ACC PRESENT(bc_h, w_p) 746 DO m = 1, bc_h(l)%ns 747 i = bc_h(l)%i(m) 748 j = bc_h(l)%j(m) 749 k = bc_h(l)%k(m) 750 w_p(k+bc_h(l)%koff,j,i) = 0.0_wp 751 ENDDO 752 ENDDO 753 754 ! 755 !-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings. 756 IF ( ibc_uv_t == 0 ) THEN 757 !$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init) 758 u_p(nzt+1,:,:) = u_init(nzt+1) 759 v_p(nzt+1,:,:) = v_init(nzt+1) 760 !$ACC END KERNELS 761 ELSEIF ( ibc_uv_t == 1 ) THEN 762 u_p(nzt+1,:,:) = u_p(nzt,:,:) 763 v_p(nzt+1,:,:) = v_p(nzt,:,:) 764 ENDIF 765 766 ! 767 !-- Vertical nesting: Vertical velocity not zero at the top of the fine grid 768 IF ( .NOT. child_domain .AND. .NOT. nesting_offline .AND. & 769 TRIM(coupling_mode) /= 'vnested_fine' ) THEN 770 !$ACC KERNELS PRESENT(w_p) 771 w_p(nzt:nzt+1,:,:) = 0.0_wp !< nzt is not a prognostic level (but cf. pres) 772 !$ACC END KERNELS 773 ENDIF 774 775 ! 776 !-- Temperature at bottom and top boundary. 777 !-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by 778 !-- the sea surface temperature of the coupled ocean model. 779 !-- Dirichlet 780 IF ( .NOT. neutral ) THEN 781 IF ( ibc_pt_b == 0 ) THEN 782 DO l = 0, 1 783 !$OMP PARALLEL DO PRIVATE( i, j, k ) 784 DO m = 1, bc_h(l)%ns 785 i = bc_h(l)%i(m) 786 j = bc_h(l)%j(m) 787 k = bc_h(l)%k(m) 788 pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i) 789 ENDDO 790 ENDDO 791 ! 792 !-- Neumann, zero-gradient 793 ELSEIF ( ibc_pt_b == 1 ) THEN 794 DO l = 0, 1 795 !$OMP PARALLEL DO PRIVATE( i, j, k ) 796 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 797 !$ACC PRESENT(bc_h, pt_p) 798 DO m = 1, bc_h(l)%ns 799 i = bc_h(l)%i(m) 800 j = bc_h(l)%j(m) 801 k = bc_h(l)%k(m) 802 pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i) 803 ENDDO 804 ENDDO 805 ENDIF 806 807 ! 808 !-- Temperature at top boundary 809 IF ( ibc_pt_t == 0 ) THEN 810 pt_p(nzt+1,:,:) = pt(nzt+1,:,:) 811 ! 812 !-- In case of nudging adjust top boundary to pt which is 813 !-- read in from NUDGING-DATA 814 IF ( nudging ) THEN 815 pt_p(nzt+1,:,:) = pt_init(nzt+1) 816 ENDIF 817 ELSEIF ( ibc_pt_t == 1 ) THEN 818 pt_p(nzt+1,:,:) = pt_p(nzt,:,:) 819 ELSEIF ( ibc_pt_t == 2 ) THEN 820 !$ACC KERNELS PRESENT(pt_p, dzu) 821 pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1) 822 !$ACC END KERNELS 823 ENDIF 824 ENDIF 825 ! 826 !-- Boundary conditions for total water content, 827 !-- bottom and top boundary (see also temperature) 828 IF ( humidity ) THEN 829 ! 830 !-- Surface conditions for constant_humidity_flux 831 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype 832 !-- the k coordinate belongs to the atmospheric grid point, therefore, set 833 !-- q_p at k-1 834 IF ( ibc_q_b == 0 ) THEN 835 836 DO l = 0, 1 837 !$OMP PARALLEL DO PRIVATE( i, j, k ) 838 DO m = 1, bc_h(l)%ns 839 i = bc_h(l)%i(m) 840 j = bc_h(l)%j(m) 841 k = bc_h(l)%k(m) 842 q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i) 843 ENDDO 844 ENDDO 845 846 ELSE 847 848 DO l = 0, 1 849 !$OMP PARALLEL DO PRIVATE( i, j, k ) 850 DO m = 1, bc_h(l)%ns 851 i = bc_h(l)%i(m) 852 j = bc_h(l)%j(m) 853 k = bc_h(l)%k(m) 854 q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i) 855 ENDDO 856 ENDDO 857 ENDIF 858 ! 859 !-- Top boundary 860 IF ( ibc_q_t == 0 ) THEN 861 q_p(nzt+1,:,:) = q(nzt+1,:,:) 862 ELSEIF ( ibc_q_t == 1 ) THEN 863 q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1) 864 ENDIF 865 ENDIF 866 ! 867 !-- Boundary conditions for scalar, 868 !-- bottom and top boundary (see also temperature) 869 IF ( passive_scalar ) THEN 870 ! 871 !-- Surface conditions for constant_humidity_flux 872 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype 873 !-- the k coordinate belongs to the atmospheric grid point, therefore, set 874 !-- s_p at k-1 875 IF ( ibc_s_b == 0 ) THEN 876 877 DO l = 0, 1 878 !$OMP PARALLEL DO PRIVATE( i, j, k ) 879 DO m = 1, bc_h(l)%ns 880 i = bc_h(l)%i(m) 881 j = bc_h(l)%j(m) 882 k = bc_h(l)%k(m) 883 s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i) 884 ENDDO 885 ENDDO 886 887 ELSE 888 889 DO l = 0, 1 890 !$OMP PARALLEL DO PRIVATE( i, j, k ) 891 DO m = 1, bc_h(l)%ns 892 i = bc_h(l)%i(m) 893 j = bc_h(l)%j(m) 894 k = bc_h(l)%k(m) 895 s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i) 896 ENDDO 897 ENDDO 898 ENDIF 899 ! 900 !-- Top boundary condition 901 IF ( ibc_s_t == 0 ) THEN 902 s_p(nzt+1,:,:) = s(nzt+1,:,:) 903 ELSEIF ( ibc_s_t == 1 ) THEN 904 s_p(nzt+1,:,:) = s_p(nzt,:,:) 905 ELSEIF ( ibc_s_t == 2 ) THEN 906 s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1) 907 ENDIF 908 909 ENDIF 910 ! 911 !-- In case of inflow or nest boundary at the south boundary the boundary for v 912 !-- is at nys and in case of inflow or nest boundary at the left boundary the 913 !-- boundary for u is at nxl. Since in prognostic_equations (cache optimized 914 !-- version) these levels are handled as a prognostic level, boundary values 915 !-- have to be restored here. 916 IF ( bc_dirichlet_s ) THEN 917 v_p(:,nys,:) = v_p(:,nys-1,:) 918 ELSEIF ( bc_dirichlet_l ) THEN 919 u_p(:,:,nxl) = u_p(:,:,nxl-1) 920 ENDIF 921 922 ! 923 !-- The same restoration for u at i=nxl and v at j=nys as above must be made 924 !-- in case of nest boundaries. This must not be done in case of vertical nesting 925 !-- mode as in that case the lateral boundaries are actually cyclic. 926 !-- Lateral oundary conditions for TKE and dissipation are set 927 !-- in tcm_boundary_conds. 928 IF ( nesting_mode /= 'vertical' .OR. nesting_offline ) THEN 929 IF ( bc_dirichlet_s ) THEN 930 v_p(:,nys,:) = v_p(:,nys-1,:) 931 ENDIF 932 IF ( bc_dirichlet_l ) THEN 933 u_p(:,:,nxl) = u_p(:,:,nxl-1) 934 ENDIF 935 ENDIF 936 937 ! 938 !-- Lateral boundary conditions for scalar quantities at the outflow. 939 !-- Lateral oundary conditions for TKE and dissipation are set 940 !-- in tcm_boundary_conds. 941 IF ( bc_radiation_s ) THEN 942 pt_p(:,nys-1,:) = pt_p(:,nys,:) 943 IF ( humidity ) THEN 944 q_p(:,nys-1,:) = q_p(:,nys,:) 945 ENDIF 946 IF ( passive_scalar ) s_p(:,nys-1,:) = s_p(:,nys,:) 947 ELSEIF ( bc_radiation_n ) THEN 948 pt_p(:,nyn+1,:) = pt_p(:,nyn,:) 949 IF ( humidity ) THEN 950 q_p(:,nyn+1,:) = q_p(:,nyn,:) 951 ENDIF 952 IF ( passive_scalar ) s_p(:,nyn+1,:) = s_p(:,nyn,:) 953 ELSEIF ( bc_radiation_l ) THEN 954 pt_p(:,:,nxl-1) = pt_p(:,:,nxl) 955 IF ( humidity ) THEN 956 q_p(:,:,nxl-1) = q_p(:,:,nxl) 957 ENDIF 958 IF ( passive_scalar ) s_p(:,:,nxl-1) = s_p(:,:,nxl) 959 ELSEIF ( bc_radiation_r ) THEN 960 pt_p(:,:,nxr+1) = pt_p(:,:,nxr) 961 IF ( humidity ) THEN 962 q_p(:,:,nxr+1) = q_p(:,:,nxr) 963 ENDIF 964 IF ( passive_scalar ) s_p(:,:,nxr+1) = s_p(:,:,nxr) 965 ENDIF 966 967 ! 968 !-- Radiation boundary conditions for the velocities at the respective outflow. 969 !-- The phase velocity is either assumed to the maximum phase velocity that 970 !-- ensures numerical stability (CFL-condition) or calculated after 971 !-- Orlanski(1976) and averaged along the outflow boundary. 972 IF ( bc_radiation_s ) THEN 973 974 IF ( use_cmax ) THEN 975 u_p(:,-1,:) = u(:,0,:) 976 v_p(:,0,:) = v(:,1,:) 977 w_p(:,-1,:) = w(:,0,:) 978 ELSEIF ( .NOT. use_cmax ) THEN 979 980 c_max = dy / dt_3d 981 982 c_u_m_l = 0.0_wp 983 c_v_m_l = 0.0_wp 984 c_w_m_l = 0.0_wp 985 986 c_u_m = 0.0_wp 987 c_v_m = 0.0_wp 988 c_w_m = 0.0_wp 989 990 ! 991 !-- Calculate the phase speeds for u, v, and w, first local and then 992 !-- average along the outflow boundary. 993 DO k = nzb+1, nzt+1 994 DO i = nxl, nxr 995 996 denom = u_m_s(k,0,i) - u_m_s(k,1,i) 997 998 IF ( denom /= 0.0_wp ) THEN 999 c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) ) 1000 IF ( c_u(k,i) < 0.0_wp ) THEN 1001 c_u(k,i) = 0.0_wp 1002 ELSEIF ( c_u(k,i) > c_max ) THEN 1003 c_u(k,i) = c_max 1004 ENDIF 1005 ELSE 1006 c_u(k,i) = c_max 1007 ENDIF 1008 1009 denom = v_m_s(k,1,i) - v_m_s(k,2,i) 1010 1011 IF ( denom /= 0.0_wp ) THEN 1012 c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) ) 1013 IF ( c_v(k,i) < 0.0_wp ) THEN 1014 c_v(k,i) = 0.0_wp 1015 ELSEIF ( c_v(k,i) > c_max ) THEN 1016 c_v(k,i) = c_max 1017 ENDIF 1018 ELSE 1019 c_v(k,i) = c_max 1020 ENDIF 1021 1022 denom = w_m_s(k,0,i) - w_m_s(k,1,i) 1023 1024 IF ( denom /= 0.0_wp ) THEN 1025 c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) ) 1026 IF ( c_w(k,i) < 0.0_wp ) THEN 1027 c_w(k,i) = 0.0_wp 1028 ELSEIF ( c_w(k,i) > c_max ) THEN 1029 c_w(k,i) = c_max 1030 ENDIF 1031 ELSE 1032 c_w(k,i) = c_max 1033 ENDIF 1034 1035 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) 1036 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) 1037 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) 1038 1039 ENDDO 1040 ENDDO 1041 1042 #if defined( __parallel ) 1043 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1044 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 1045 MPI_SUM, comm1dx, ierr ) 1046 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1047 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 1048 MPI_SUM, comm1dx, ierr ) 1049 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1050 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 1051 MPI_SUM, comm1dx, ierr ) 1052 #else 1053 c_u_m = c_u_m_l 1054 c_v_m = c_v_m_l 1055 c_w_m = c_w_m_l 1056 #endif 1057 1058 c_u_m = c_u_m / (nx+1) 1059 c_v_m = c_v_m / (nx+1) 1060 c_w_m = c_w_m / (nx+1) 1061 1062 ! 1063 !-- Save old timelevels for the next timestep 1064 IF ( intermediate_timestep_count == 1 ) THEN 1065 u_m_s(:,:,:) = u(:,0:1,:) 1066 v_m_s(:,:,:) = v(:,1:2,:) 1067 w_m_s(:,:,:) = w(:,0:1,:) 1068 ENDIF 1069 1070 ! 1071 !-- Calculate the new velocities 1072 DO k = nzb+1, nzt+1 1073 DO i = nxlg, nxrg 1074 u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * & 1075 ( u(k,-1,i) - u(k,0,i) ) * ddy 1076 1077 v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * & 1078 ( v(k,0,i) - v(k,1,i) ) * ddy 1079 1080 w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * & 1081 ( w(k,-1,i) - w(k,0,i) ) * ddy 1082 ENDDO 1083 ENDDO 1084 1085 ! 1086 !-- Bottom boundary at the outflow 1087 IF ( ibc_uv_b == 0 ) THEN 1088 u_p(nzb,-1,:) = 0.0_wp 1089 v_p(nzb,0,:) = 0.0_wp 1090 ELSE 1091 u_p(nzb,-1,:) = u_p(nzb+1,-1,:) 1092 v_p(nzb,0,:) = v_p(nzb+1,0,:) 1093 ENDIF 1094 w_p(nzb,-1,:) = 0.0_wp 1095 1096 ! 1097 !-- Top boundary at the outflow 1098 IF ( ibc_uv_t == 0 ) THEN 1099 u_p(nzt+1,-1,:) = u_init(nzt+1) 1100 v_p(nzt+1,0,:) = v_init(nzt+1) 1101 ELSE 1102 u_p(nzt+1,-1,:) = u_p(nzt,-1,:) 1103 v_p(nzt+1,0,:) = v_p(nzt,0,:) 1104 ENDIF 1105 w_p(nzt:nzt+1,-1,:) = 0.0_wp 1106 1107 ENDIF 1108 1109 ENDIF 1110 1111 IF ( bc_radiation_n ) THEN 1112 1113 IF ( use_cmax ) THEN 1114 u_p(:,ny+1,:) = u(:,ny,:) 1115 v_p(:,ny+1,:) = v(:,ny,:) 1116 w_p(:,ny+1,:) = w(:,ny,:) 1117 ELSEIF ( .NOT. use_cmax ) THEN 1118 1119 c_max = dy / dt_3d 1120 1121 c_u_m_l = 0.0_wp 1122 c_v_m_l = 0.0_wp 1123 c_w_m_l = 0.0_wp 1124 1125 c_u_m = 0.0_wp 1126 c_v_m = 0.0_wp 1127 c_w_m = 0.0_wp 1128 1129 ! 1130 !-- Calculate the phase speeds for u, v, and w, first local and then 1131 !-- average along the outflow boundary. 1132 DO k = nzb+1, nzt+1 1133 DO i = nxl, nxr 1134 1135 denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i) 1136 1137 IF ( denom /= 0.0_wp ) THEN 1138 c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) ) 1139 IF ( c_u(k,i) < 0.0_wp ) THEN 1140 c_u(k,i) = 0.0_wp 1141 ELSEIF ( c_u(k,i) > c_max ) THEN 1142 c_u(k,i) = c_max 1143 ENDIF 1144 ELSE 1145 c_u(k,i) = c_max 1146 ENDIF 1147 1148 denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i) 1149 1150 IF ( denom /= 0.0_wp ) THEN 1151 c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) ) 1152 IF ( c_v(k,i) < 0.0_wp ) THEN 1153 c_v(k,i) = 0.0_wp 1154 ELSEIF ( c_v(k,i) > c_max ) THEN 1155 c_v(k,i) = c_max 1156 ENDIF 1157 ELSE 1158 c_v(k,i) = c_max 1159 ENDIF 1160 1161 denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i) 1162 1163 IF ( denom /= 0.0_wp ) THEN 1164 c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) ) 1165 IF ( c_w(k,i) < 0.0_wp ) THEN 1166 c_w(k,i) = 0.0_wp 1167 ELSEIF ( c_w(k,i) > c_max ) THEN 1168 c_w(k,i) = c_max 1169 ENDIF 1170 ELSE 1171 c_w(k,i) = c_max 1172 ENDIF 1173 1174 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) 1175 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) 1176 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) 1177 1178 ENDDO 1179 ENDDO 1180 1181 #if defined( __parallel ) 1182 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1183 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 1184 MPI_SUM, comm1dx, ierr ) 1185 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1186 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 1187 MPI_SUM, comm1dx, ierr ) 1188 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 1189 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 1190 MPI_SUM, comm1dx, ierr ) 1191 #else 1192 c_u_m = c_u_m_l 1193 c_v_m = c_v_m_l 1194 c_w_m = c_w_m_l 1195 #endif 1196 1197 c_u_m = c_u_m / (nx+1) 1198 c_v_m = c_v_m / (nx+1) 1199 c_w_m = c_w_m / (nx+1) 1200 1201 ! 1202 !-- Save old timelevels for the next timestep 1203 IF ( intermediate_timestep_count == 1 ) THEN 1204 u_m_n(:,:,:) = u(:,ny-1:ny,:) 1205 v_m_n(:,:,:) = v(:,ny-1:ny,:) 1206 w_m_n(:,:,:) = w(:,ny-1:ny,:) 1207 ENDIF 1208 1209 ! 1210 !-- Calculate the new velocities 1211 DO k = nzb+1, nzt+1 1212 DO i = nxlg, nxrg 1213 u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * & 1214 ( u(k,ny+1,i) - u(k,ny,i) ) * ddy 1215 1216 v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * & 1217 ( v(k,ny+1,i) - v(k,ny,i) ) * ddy 1218 1219 w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * & 1220 ( w(k,ny+1,i) - w(k,ny,i) ) * ddy 1221 ENDDO 1222 ENDDO 1223 1224 ! 1225 !-- Bottom boundary at the outflow 1226 IF ( ibc_uv_b == 0 ) THEN 1227 u_p(nzb,ny+1,:) = 0.0_wp 1228 v_p(nzb,ny+1,:) = 0.0_wp 1229 ELSE 1230 u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) 1231 v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:) 1232 ENDIF 1233 w_p(nzb,ny+1,:) = 0.0_wp 1234 1235 ! 1236 !-- Top boundary at the outflow 1237 IF ( ibc_uv_t == 0 ) THEN 1238 u_p(nzt+1,ny+1,:) = u_init(nzt+1) 1239 v_p(nzt+1,ny+1,:) = v_init(nzt+1) 1240 ELSE 1241 u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:) 1242 v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:) 1243 ENDIF 1244 w_p(nzt:nzt+1,ny+1,:) = 0.0_wp 1245 1246 ENDIF 1247 1248 ENDIF 1249 1250 IF ( bc_radiation_l ) THEN 1251 1252 IF ( use_cmax ) THEN 1253 u_p(:,:,0) = u(:,:,1) 1254 v_p(:,:,-1) = v(:,:,0) 1255 w_p(:,:,-1) = w(:,:,0) 1256 ELSEIF ( .NOT. use_cmax ) THEN 1257 1258 c_max = dx / dt_3d 1259 1260 c_u_m_l = 0.0_wp 1261 c_v_m_l = 0.0_wp 1262 c_w_m_l = 0.0_wp 1263 1264 c_u_m = 0.0_wp 1265 c_v_m = 0.0_wp 1266 c_w_m = 0.0_wp 1267 1268 ! 1269 !-- Calculate the phase speeds for u, v, and w, first local and then 1270 !-- average along the outflow boundary. 1271 DO k = nzb+1, nzt+1 1272 DO j = nys, nyn 1273 1274 denom = u_m_l(k,j,1) - u_m_l(k,j,2) 1275 1276 IF ( denom /= 0.0_wp ) THEN 1277 c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) ) 1278 IF ( c_u(k,j) < 0.0_wp ) THEN 1279 c_u(k,j) = 0.0_wp 1280 ELSEIF ( c_u(k,j) > c_max ) THEN 1281 c_u(k,j) = c_max 1282 ENDIF 1283 ELSE 1284 c_u(k,j) = c_max 1285 ENDIF 1286 1287 denom = v_m_l(k,j,0) - v_m_l(k,j,1) 1288 1289 IF ( denom /= 0.0_wp ) THEN 1290 c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) ) 1291 IF ( c_v(k,j) < 0.0_wp ) THEN 1292 c_v(k,j) = 0.0_wp 1293 ELSEIF ( c_v(k,j) > c_max ) THEN 1294 c_v(k,j) = c_max 1295 ENDIF 1296 ELSE 1297 c_v(k,j) = c_max 1298 ENDIF 1299 1300 denom = w_m_l(k,j,0) - w_m_l(k,j,1) 1301 1302 IF ( denom /= 0.0_wp ) THEN 1303 c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) ) 1304 IF ( c_w(k,j) < 0.0_wp ) THEN 1305 c_w(k,j) = 0.0_wp 1306 ELSEIF ( c_w(k,j) > c_max ) THEN 1307 c_w(k,j) = c_max 1308 ENDIF 1309 ELSE 1310 c_w(k,j) = c_max 1311 ENDIF 1312 1313 c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) 1314 c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) 1315 c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) 1316 1317 ENDDO 1318 ENDDO 1319 1320 #if defined( __parallel ) 1321 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1322 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 1323 MPI_SUM, comm1dy, ierr ) 1324 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1325 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 1326 MPI_SUM, comm1dy, ierr ) 1327 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1328 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 1329 MPI_SUM, comm1dy, ierr ) 1330 #else 1331 c_u_m = c_u_m_l 1332 c_v_m = c_v_m_l 1333 c_w_m = c_w_m_l 1334 #endif 1335 1336 c_u_m = c_u_m / (ny+1) 1337 c_v_m = c_v_m / (ny+1) 1338 c_w_m = c_w_m / (ny+1) 1339 1340 ! 1341 !-- Save old timelevels for the next timestep 1342 IF ( intermediate_timestep_count == 1 ) THEN 1343 u_m_l(:,:,:) = u(:,:,1:2) 1344 v_m_l(:,:,:) = v(:,:,0:1) 1345 w_m_l(:,:,:) = w(:,:,0:1) 1346 ENDIF 1347 1348 ! 1349 !-- Calculate the new velocities 1350 DO k = nzb+1, nzt+1 1351 DO j = nysg, nyng 1352 u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * & 1353 ( u(k,j,0) - u(k,j,1) ) * ddx 1354 1355 v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * & 1356 ( v(k,j,-1) - v(k,j,0) ) * ddx 1357 1358 w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * & 1359 ( w(k,j,-1) - w(k,j,0) ) * ddx 1360 ENDDO 1361 ENDDO 1362 1363 ! 1364 !-- Bottom boundary at the outflow 1365 IF ( ibc_uv_b == 0 ) THEN 1366 u_p(nzb,:,0) = 0.0_wp 1367 v_p(nzb,:,-1) = 0.0_wp 1368 ELSE 1369 u_p(nzb,:,0) = u_p(nzb+1,:,0) 1370 v_p(nzb,:,-1) = v_p(nzb+1,:,-1) 1371 ENDIF 1372 w_p(nzb,:,-1) = 0.0_wp 1373 1374 ! 1375 !-- Top boundary at the outflow 1376 IF ( ibc_uv_t == 0 ) THEN 1377 u_p(nzt+1,:,0) = u_init(nzt+1) 1378 v_p(nzt+1,:,-1) = v_init(nzt+1) 1379 ELSE 1380 u_p(nzt+1,:,0) = u_p(nzt,:,0) 1381 v_p(nzt+1,:,-1) = v_p(nzt,:,-1) 1382 ENDIF 1383 w_p(nzt:nzt+1,:,-1) = 0.0_wp 1384 1385 ENDIF 1386 1387 ENDIF 1388 1389 IF ( bc_radiation_r ) THEN 1390 1391 IF ( use_cmax ) THEN 1392 u_p(:,:,nx+1) = u(:,:,nx) 1393 v_p(:,:,nx+1) = v(:,:,nx) 1394 w_p(:,:,nx+1) = w(:,:,nx) 1395 ELSEIF ( .NOT. use_cmax ) THEN 1396 1397 c_max = dx / dt_3d 1398 1399 c_u_m_l = 0.0_wp 1400 c_v_m_l = 0.0_wp 1401 c_w_m_l = 0.0_wp 1402 1403 c_u_m = 0.0_wp 1404 c_v_m = 0.0_wp 1405 c_w_m = 0.0_wp 1406 1407 ! 1408 !-- Calculate the phase speeds for u, v, and w, first local and then 1409 !-- average along the outflow boundary. 1410 DO k = nzb+1, nzt+1 1411 DO j = nys, nyn 1412 1413 denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1) 1414 1415 IF ( denom /= 0.0_wp ) THEN 1416 c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) ) 1417 IF ( c_u(k,j) < 0.0_wp ) THEN 1418 c_u(k,j) = 0.0_wp 1419 ELSEIF ( c_u(k,j) > c_max ) THEN 1420 c_u(k,j) = c_max 1421 ENDIF 1422 ELSE 1423 c_u(k,j) = c_max 1424 ENDIF 1425 1426 denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1) 1427 1428 IF ( denom /= 0.0_wp ) THEN 1429 c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) ) 1430 IF ( c_v(k,j) < 0.0_wp ) THEN 1431 c_v(k,j) = 0.0_wp 1432 ELSEIF ( c_v(k,j) > c_max ) THEN 1433 c_v(k,j) = c_max 1434 ENDIF 1435 ELSE 1436 c_v(k,j) = c_max 1437 ENDIF 1438 1439 denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1) 1440 1441 IF ( denom /= 0.0_wp ) THEN 1442 c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) ) 1443 IF ( c_w(k,j) < 0.0_wp ) THEN 1444 c_w(k,j) = 0.0_wp 1445 ELSEIF ( c_w(k,j) > c_max ) THEN 1446 c_w(k,j) = c_max 1447 ENDIF 1448 ELSE 1449 c_w(k,j) = c_max 1450 ENDIF 1451 1452 c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) 1453 c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) 1454 c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) 1455 1456 ENDDO 1457 ENDDO 1458 1459 #if defined( __parallel ) 1460 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1461 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 1462 MPI_SUM, comm1dy, ierr ) 1463 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1464 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 1465 MPI_SUM, comm1dy, ierr ) 1466 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 1467 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 1468 MPI_SUM, comm1dy, ierr ) 1469 #else 1470 c_u_m = c_u_m_l 1471 c_v_m = c_v_m_l 1472 c_w_m = c_w_m_l 1473 #endif 1474 1475 c_u_m = c_u_m / (ny+1) 1476 c_v_m = c_v_m / (ny+1) 1477 c_w_m = c_w_m / (ny+1) 1478 1479 ! 1480 !-- Save old timelevels for the next timestep 1481 IF ( intermediate_timestep_count == 1 ) THEN 1482 u_m_r(:,:,:) = u(:,:,nx-1:nx) 1483 v_m_r(:,:,:) = v(:,:,nx-1:nx) 1484 w_m_r(:,:,:) = w(:,:,nx-1:nx) 1485 ENDIF 1486 1487 ! 1488 !-- Calculate the new velocities 1489 DO k = nzb+1, nzt+1 1490 DO j = nysg, nyng 1491 u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * & 1492 ( u(k,j,nx+1) - u(k,j,nx) ) * ddx 1493 1494 v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * & 1495 ( v(k,j,nx+1) - v(k,j,nx) ) * ddx 1496 1497 w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * & 1498 ( w(k,j,nx+1) - w(k,j,nx) ) * ddx 1499 ENDDO 1500 ENDDO 1501 1502 ! 1503 !-- Bottom boundary at the outflow 1504 IF ( ibc_uv_b == 0 ) THEN 1505 u_p(nzb,:,nx+1) = 0.0_wp 1506 v_p(nzb,:,nx+1) = 0.0_wp 1507 ELSE 1508 u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) 1509 v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1) 1510 ENDIF 1511 w_p(nzb,:,nx+1) = 0.0_wp 1512 1513 ! 1514 !-- Top boundary at the outflow 1515 IF ( ibc_uv_t == 0 ) THEN 1516 u_p(nzt+1,:,nx+1) = u_init(nzt+1) 1517 v_p(nzt+1,:,nx+1) = v_init(nzt+1) 1518 ELSE 1519 u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1) 1520 v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1) 1521 ENDIF 1522 w_p(nzt:nzt+1,:,nx+1) = 0.0_wp 1523 1524 ENDIF 1525 1526 ENDIF 1527 1528 END SUBROUTINE dynamics_boundary_conditions 659 1529 !------------------------------------------------------------------------------! 660 1530 ! Description: -
palm/trunk/SOURCE/module_interface.f90
r4272 r4281 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Added dynamics boundary conditions 28 ! 29 ! 4272 2019-10-23 15:18:57Z schwenkel 27 30 ! Further modularization of boundary conditions: moved boundary conditions to 28 31 ! respective modules … … 179 182 dynamics_exchange_horiz, & 180 183 dynamics_prognostic_equations, & 184 dynamics_boundary_conditions, & 181 185 dynamics_swap_timelevel, & 182 186 dynamics_3d_data_averaging, & … … 197 201 tcm_actions, & 198 202 tcm_prognostic_equations, & 203 tcm_boundary_conds, & 199 204 tcm_swap_timelevel, & 200 205 tcm_3d_data_averaging, & … … 1299 1304 IF ( debug_output_timestep ) CALL debug_message( 'module-specific boundary_conditions', 'start' ) 1300 1305 1306 CALL dynamics_boundary_conditions 1307 CALL tcm_boundary_conds 1308 1301 1309 IF ( bulk_cloud_model ) CALL bcm_boundary_conditions 1302 1310 IF ( air_chemistry ) CALL chem_boundary_conditions -
palm/trunk/SOURCE/time_integration.f90
r4276 r4281 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Moved boundary conditions to module interface 28 ! 29 ! 4276 2019-10-28 16:03:29Z schwenkel 27 30 ! Further modularization of lpm code components 28 31 ! … … 730 733 ! 731 734 !-- Boundary conditions for the prognostic quantities (except of the 732 !-- velocities at the outflow in case of a non-cyclic lateral wall) 733 CALL boundary_conds 734 ! 735 !-- Boundary conditions for module-specific variables 735 !-- velocities at the outflow in case of a non-cyclic lateral wall) and 736 !-- boundary conditions for module-specific variables 736 737 CALL module_interface_boundary_conditions 737 738 !
Note: See TracChangeset
for help on using the changeset viewer.
