Changeset 4281 for palm/trunk/SOURCE/dynamics_mod.f90
 Timestamp:
 Oct 29, 2019 3:15:39 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/dynamics_mod.f90
r4097 r4281 25 25 !  26 26 ! $Id$ 27 ! Moved boundary conditions in dynamics module 28 ! 29 ! 4097 20190715 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 downwardfacing 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 downwardfacing 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, zerogradient 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 NUDGINGDATA 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 nonnatural and natural walls. Note, in walldatatype 832 ! the k coordinate belongs to the atmospheric grid point, therefore, set 833 ! q_p at k1 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 nonnatural and natural walls. Note, in walldatatype 873 ! the k coordinate belongs to the atmospheric grid point, therefore, set 874 ! s_p at k1 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(:,nys1,:) 918 ELSEIF ( bc_dirichlet_l ) THEN 919 u_p(:,:,nxl) = u_p(:,:,nxl1) 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(:,nys1,:) 931 ENDIF 932 IF ( bc_dirichlet_l ) THEN 933 u_p(:,:,nxl) = u_p(:,:,nxl1) 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(:,nys1,:) = pt_p(:,nys,:) 943 IF ( humidity ) THEN 944 q_p(:,nys1,:) = q_p(:,nys,:) 945 ENDIF 946 IF ( passive_scalar ) s_p(:,nys1,:) = 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(:,:,nxl1) = pt_p(:,:,nxl) 955 IF ( humidity ) THEN 956 q_p(:,:,nxl1) = q_p(:,:,nxl) 957 ENDIF 958 IF ( passive_scalar ) s_p(:,:,nxl1) = 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 (CFLcondition) 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), nztnzb, 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), nztnzb, 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), nztnzb, 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,ny1,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,ny1,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,ny1,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), nztnzb, 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), nztnzb, 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), nztnzb, 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(:,ny1:ny,:) 1205 v_m_n(:,:,:) = v(:,ny1:ny,:) 1206 w_m_n(:,:,:) = w(:,ny1: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), nztnzb, 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), nztnzb, 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), nztnzb, 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,nx1) 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,nx1) 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,nx1) 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), nztnzb, 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), nztnzb, 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), nztnzb, 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(:,:,nx1:nx) 1483 v_m_r(:,:,:) = v(:,:,nx1:nx) 1484 w_m_r(:,:,:) = w(:,:,nx1: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:
Note: See TracChangeset
for help on using the changeset viewer.