Changeset 3708 for palm/trunk
 Timestamp:
 Jan 30, 2019 12:58:13 PM (6 years ago)
 Location:
 palm/trunk
 Files:

 8 deleted
 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/pmc_interface_mod.f90
r3697 r3708 25 25 !  26 26 ! $Id$ 27 ! Checks for parent / child grid line matching introduced. 28 ! Interpolation of nestboundarytangential velocity components revised. 29 ! 30 ! 3697 20190124 17:16:13Z hellstea 27 31 ! Bugfix: upper kbound in the child initialization interpolation 28 32 ! pmci_interp_1sto_all corrected. … … 1475 1479 ! included in the interpolation algorithms. 1476 1480 CALL pmci_init_anterp_tophat 1481 ! 1482 ! Check that the child and parent grid lines match 1483 CALL pmci_check_grid_matching 1477 1484 1478 1485 ENDIF … … 1928 1935 CEILING( parentdzmax / dzmin ) 1929 1936 ALLOCATE( celltmpd(1:acsize) ) 1930 1937 1931 1938 END SUBROUTINE pmci_init_interp_tril 1932 1939 … … 1941 1948 igsr = NINT( cg%dx / dx, iwp ) 1942 1949 jgsr = NINT( cg%dy / dy, iwp ) 1943 kgsr = NINT( cg%dz u(cg%nz+1) / dzu(nzt+1), iwp )1950 kgsr = NINT( cg%dzw(1) / dzw(1), iwp ) 1944 1951 write(9,"('igsr, jgsr, kgsr: ',3(i3,2x))") igsr, jgsr, kgsr 1945 1952 flush(9) … … 3899 3906 END SUBROUTINE pmci_init_anterp_tophat 3900 3907 3901 #endif 3908 3909 3910 SUBROUTINE pmci_check_grid_matching 3911 ! 3912 ! Check that the grid lines of child and parent do match 3913 IMPLICIT NONE 3914 REAL(wp), PARAMETER :: tolefac = 1.0E6_wp 3915 REAL(wp) :: tolex 3916 REAL(wp) :: toley 3917 REAL(wp) :: tolez 3918 INTEGER(iwp) :: non_matching_lower_left_corner 3919 INTEGER(iwp) :: non_int_gsr_x 3920 INTEGER(iwp) :: non_int_gsr_y 3921 INTEGER(iwp) :: non_int_gsr_z 3922 3923 3924 non_matching_lower_left_corner = 0 3925 non_int_gsr_x = 0 3926 non_int_gsr_y = 0 3927 non_int_gsr_z = 0 3928 3929 IF ( myid == 0 ) THEN 3930 3931 tolex = tolefac * dx 3932 toley = tolefac * dy 3933 tolez = tolefac * MINVAL( dzw ) 3934 ! 3935 ! First check that the child grid lower left corner matches the paren grid lines. 3936 IF ( MOD( lower_left_coord_x, cg%dx ) > tolex ) non_matching_lower_left_corner = 1 3937 IF ( MOD( lower_left_coord_y, cg%dy ) > toley ) non_matching_lower_left_corner = 1 3938 ! 3939 ! Then check that the gridspacing ratios in each direction are integer valued. 3940 IF ( MOD( cg%dx, dx ) > tolex ) non_int_gsr_x = 1 3941 IF ( MOD( cg%dy, dy ) > toley ) non_int_gsr_y = 1 3942 DO n = 0, kctw+1 3943 IF ( ABS( cg%zw(n)  zw(kflw(n)) ) > tolez ) non_int_gsr_z = 1 3944 ENDDO 3945 3946 ENDIF 3947 3948 CALL MPI_BCAST( non_matching_lower_left_corner, 1, MPI_INTEGER, 0, & 3949 comm2d, ierr ) 3950 IF ( non_matching_lower_left_corner > 0 ) THEN 3951 WRITE ( message_string, * ) 'nested child domain lower left ', & 3952 'corner must match its parent grid lines' 3953 CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 ) 3954 ENDIF 3955 3956 CALL MPI_BCAST( non_int_gsr_x, 1, MPI_INTEGER, 0, comm2d, ierr ) 3957 IF ( non_int_gsr_x > 0 ) THEN 3958 WRITE ( message_string, * ) 'nesting gridspacing ratio ', & 3959 '( parent dx / child dx ) must have an integer value' 3960 CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 ) 3961 ENDIF 3962 3963 CALL MPI_BCAST( non_int_gsr_y, 1, MPI_INTEGER, 0, comm2d, ierr ) 3964 IF ( non_int_gsr_y > 0 ) THEN 3965 WRITE ( message_string, * ) 'nesting gridspacing ratio ', & 3966 '( parent dy / child dy ) must have an integer value' 3967 CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 ) 3968 ENDIF 3969 3970 CALL MPI_BCAST( non_int_gsr_z, 1, MPI_INTEGER, 0, comm2d, ierr ) 3971 IF ( non_int_gsr_z > 0 ) THEN 3972 WRITE ( message_string, * ) 'nesting gridspacing ratio ', & 3973 '( parent dz / child dz ) must have an integer value for each zlevel' 3974 CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 ) 3975 ENDIF 3976 3977 END SUBROUTINE pmci_check_grid_matching 3978 3979 3980 #endif 3902 3981 END SUBROUTINE pmci_setup_child 3903 3982 … … 4557 4636 DO m = mb, me 4558 4637 DO n = 0, kct + 1 4559 4638 4560 4639 DO i = ifl(l), ifu(l) 4561 4640 DO j = jfl(m), jfl(m+1)1 … … 4594 4673 DO m = mb, me 4595 4674 DO n = 0, kct + 1 4596 4675 4597 4676 DO i = ifl(l), ifu(l) 4598 DO j = jfl(m), jfu(m) 4677 DO j = jfl(m), jfu(m) 4599 4678 DO k = kfl(n), MIN( kfu(n), nzt+1 ) 4600 4679 f(k,j,i) = fc(n,m,l) … … 5785 5864 REAL(wp) :: rcorr !< Average reversibility correction for the whole anterpolation cell 5786 5865 REAL(wp) :: rcorr_ijk !< Reversibility correction distributed to the individual childgrid nodes 5866 5867 real(wp), parameter :: c1 = 2.0_wp / 6.0_wp 5868 real(wp), parameter :: c2 = 5.0_wp / 6.0_wp 5869 real(wp), parameter :: c3 = 1.0_wp / 6.0_wp 5787 5870 5788 5871 ! … … 5858 5941 ENDIF 5859 5942 5860 IF ( var == 'u' ) THEN 5861 !AH! 5862 !AH! Substitute the necessary parentgrid data to the work array workarrc_lr. 5863 !AH workarrc_lr = 0.0_wp 5864 !AH IF ( pdims(2) > 1 ) THEN 5865 !AH#if defined( __parallel ) 5866 !AH IF ( nys == 0 ) THEN 5867 !AH workarrc_lr(0:cg%nz+1,jcsw:jcnw1,0:2) & 5868 !AH = fc(0:cg%nz+1,jcsw:jcnw1,lbeg:lbeg+2) 5869 !AH ELSE IF ( nyn == ny ) THEN 5870 !AH workarrc_lr(0:cg%nz+1,jcsw+1:jcnw,0:2) & 5871 !AH = fc(0:cg%nz+1,jcsw+1:jcnw,lbeg:lbeg+2) 5872 !AH ELSE 5873 !AH workarrc_lr(0:cg%nz+1,jcsw+1:jcnw1,0:2) & 5874 !AH = fc(0:cg%nz+1,jcsw+1:jcnw1,lbeg:lbeg+2) 5875 !AH ENDIF 5876 !AH! 5877 !AH! Southnorth exchange if more than one PE subdomain in the ydirection. 5878 !AH! Note that in case of 3D nesting the south (psouth == MPI_PROC_NULL) 5879 !AH! and north (pnorth == MPI_PROC_NULL) boundaries are not exchanged 5880 !AH! because the nest domain is not cyclic. 5881 !AH! From south to north 5882 !AH CALL MPI_SENDRECV( workarrc_lr(0,jcsw+1,0), 1, & 5883 !AH workarrc_lr_exchange_type, psouth, 0, & 5884 !AH workarrc_lr(0,jcnw,0), 1, & 5885 !AH workarrc_lr_exchange_type, pnorth, 0, & 5886 !AH comm2d, status, ierr ) 5887 !AH! 5888 !AH! From north to south 5889 !AH CALL MPI_SENDRECV( workarrc_lr(0,jcnw1,0), 1, & 5890 !AH workarrc_lr_exchange_type, pnorth, 1, & 5891 !AH workarrc_lr(0,jcsw,0), 1, & 5892 !AH workarrc_lr_exchange_type, psouth, 1, & 5893 !AH comm2d, status, ierr ) 5894 !AH#endif 5895 !AH ELSE 5896 !AH workarrc_lr(0:cg%nz+1,jcsw:jcnw,0:2) & 5897 !AH = fc(0:cg%nz+1,jcsw:jcnw,lbeg:lbeg+2) 5898 !AH ENDIF 5899 5900 DO m = jcsw+1, jcnw1 5901 DO n = 0, kct 5902 5903 DO j = jfl(m), jfu(m) 5904 DO k = kfl(n), kfu(n) 5905 f(k,j,ibc) = fc(n,m,l) 5943 IF ( var == 'v' .OR. var == 'w' ) THEN 5944 ! 5945 ! Substitute the necessary parentgrid data to the work array workarrc_lr. 5946 workarrc_lr = 0.0_wp 5947 IF ( pdims(2) > 1 ) THEN 5948 #if defined( __parallel ) 5949 IF ( nys == 0 ) THEN 5950 workarrc_lr(0:cg%nz+1,jcsw:jcnw1,0:2) & 5951 = fc(0:cg%nz+1,jcsw:jcnw1,lbeg:lbeg+2) 5952 ELSE IF ( nyn == ny ) THEN 5953 workarrc_lr(0:cg%nz+1,jcsw+1:jcnw,0:2) & 5954 = fc(0:cg%nz+1,jcsw+1:jcnw,lbeg:lbeg+2) 5955 ELSE 5956 workarrc_lr(0:cg%nz+1,jcsw+1:jcnw1,0:2) & 5957 = fc(0:cg%nz+1,jcsw+1:jcnw1,lbeg:lbeg+2) 5958 ENDIF 5959 ! 5960 ! Southnorth exchange if more than one PE subdomain in the ydirection. 5961 ! Note that in case of 3D nesting the south (psouth == MPI_PROC_NULL) 5962 ! and north (pnorth == MPI_PROC_NULL) boundaries are not exchanged 5963 ! because the nest domain is not cyclic. 5964 ! From south to north 5965 CALL MPI_SENDRECV( workarrc_lr(0,jcsw+1,0), 1, & 5966 workarrc_lr_exchange_type, psouth, 0, & 5967 workarrc_lr(0,jcnw,0), 1, & 5968 workarrc_lr_exchange_type, pnorth, 0, & 5969 comm2d, status, ierr ) 5970 ! 5971 ! From north to south 5972 CALL MPI_SENDRECV( workarrc_lr(0,jcnw1,0), 1, & 5973 workarrc_lr_exchange_type, pnorth, 1, & 5974 workarrc_lr(0,jcsw,0), 1, & 5975 workarrc_lr_exchange_type, psouth, 1, & 5976 comm2d, status, ierr ) 5977 #endif 5978 ELSE 5979 workarrc_lr(0:cg%nz+1,jcsw:jcnw,0:2) & 5980 = fc(0:cg%nz+1,jcsw:jcnw,lbeg:lbeg+2) 5981 ENDIF 5982 5983 IF ( var == 'v' ) THEN 5984 5985 DO m = jcsw+1, jcnw 5986 DO n = 0, kct 5987 ! 5988 ! First substitute only the matchingnode values 5989 f(kfl(n):kfu(n),jfl(m),ibc) = workarrc_lr(n,m,lw) 5990 5991 ENDDO 5992 ENDDO 5993 5994 DO m = jcsw+1, jcnw1 5995 DO n = 0, kct 5996 ! 5997 ! Then fill up the nodes in between with the averages 5998 DO j = jfl(m)+1, jfl(m+1)1 5999 f(kfl(n):kfu(n),j,ibc) = 0.5_wp * ( f(kfl(n):kfu(n),jfl(m),ibc) & 6000 + f(kfl(n):kfu(n),jfl(m+1),ibc) ) 5906 6001 ENDDO 6002 5907 6003 ENDDO 5908 5909 6004 ENDDO 5910 ENDDO 5911 5912 ELSE IF ( var == 'v' ) THEN 5913 5914 DO m = jcsw+1, jcnw1 5915 DO n = 0, kct 5916 5917 DO j = jfl(m), jfl(m+1)1 5918 DO k = kfl(n), kfu(n) 5919 f(k,j,ibc) = fc(n,m,l) 6005 6006 ELSE IF ( var == 'w' ) THEN 6007 6008 DO m = jcsw+1, jcnw1 6009 DO n = 0, kct + 1 ! It is important to go up to kct+1 6010 ! 6011 ! First substitute only the matchingnode values 6012 f(kfu(n),jfl(m):jfu(m),ibc) = workarrc_lr(n,m,lw) 6013 6014 ENDDO 6015 ENDDO 6016 6017 DO m = jcsw+1, jcnw1 6018 DO n = 1, kct + 1 ! It is important to go up to kct+1 6019 ! 6020 ! Then fill up the nodes in between with the averages 6021 DO k = kfu(n1)+1, kfu(n)1 6022 f(k,jfl(m):jfu(m),ibc) = 0.5_wp * ( f(kfu(n1),jfl(m):jfu(m),ibc) & 6023 + f(kfu(n),jfl(m):jfu(m),ibc) ) 5920 6024 ENDDO 6025 5921 6026 ENDDO 5922 5923 6027 ENDDO 5924 ENDDO 5925 5926 ELSE IF ( var == 'w' ) THEN 5927 5928 DO m = jcsw+1, jcnw1 5929 DO n = 1, kct + 1 ! It is important to go up to kct+1 5930 5931 DO j = jfl(m), jfu(m) 5932 f(nzb,j,ibc) = 0.0_wp ! Because the nloop starts from n=1 instead of 0 5933 DO k = kfu(n1)+1, kfu(n) 5934 f(k,j,ibc) = fc(n,m,l) 5935 ENDDO 5936 ENDDO 5937 5938 ENDDO 5939 ENDDO 5940 5941 ELSE ! scalars 6028 6029 ENDIF 6030 6031 ELSE ! u or scalars 5942 6032 5943 6033 DO m = jcsw+1, jcnw1 … … 6081 6171 REAL(wp) :: rcorr !< Average reversibility correction for the whole anterpolation cell 6082 6172 REAL(wp) :: rcorr_ijk !< Reversibility correction distributed to the individual childgrid nodes 6173 6174 real(wp), parameter :: c1 = 2.0_wp / 6.0_wp 6175 real(wp), parameter :: c2 = 5.0_wp / 6.0_wp 6176 real(wp), parameter :: c3 = 1.0_wp / 6.0_wp 6083 6177 6084 6178 ! … … 6154 6248 ENDIF 6155 6249 6156 IF ( var == 'v' ) THEN 6157 !AH! 6158 !AH! Substitute the necessary parentgrid data to the work array workarrc_sn. 6159 !AH workarrc_sn = 0.0_wp 6160 !AH IF ( pdims(1) > 1 ) THEN 6161 !AH#if defined( __parallel ) 6162 !AH IF ( nxl == 0 ) THEN ! if ( bc_dirichlet_l ) 6163 !AH workarrc_sn(0:cg%nz+1,0:2,iclw:icrw1) & 6164 !AH = fc(0:cg%nz+1,mbeg:mbeg+2,iclw:icrw1) 6165 !AH ELSE IF ( nxr == nx ) THEN ! if ( bc_dirichlet_r ) 6166 !AH workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw) & 6167 !AH = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw) 6168 !AH ELSE 6169 !AH workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw1) & 6170 !AH = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw1) 6171 !AH ENDIF 6172 !AH! 6173 !AH! Leftright exchange if more than one PE subdomain in the xdirection. 6174 !AH! Note that in case of 3D nesting the left (pleft == MPI_PROC_NULL) and 6175 !AH! right (pright == MPI_PROC_NULL) boundaries are not exchanged because 6176 !AH! the nest domain is not cyclic. 6177 !AH! From left to right 6178 !AH CALL MPI_SENDRECV( workarrc_sn(0,0,iclw+1), 1, & 6179 !AH workarrc_sn_exchange_type, pleft, 0, & 6180 !AH workarrc_sn(0,0,icrw), 1, & 6181 !AH workarrc_sn_exchange_type, pright, 0, & 6182 !AH comm2d, status, ierr ) 6183 !AH! 6184 !AH! From right to left 6185 !AH CALL MPI_SENDRECV( workarrc_sn(0,0,icrw1), 1, & 6186 !AH workarrc_sn_exchange_type, pright, 1, & 6187 !AH workarrc_sn(0,0,iclw), 1, & 6188 !AH workarrc_sn_exchange_type, pleft, 1, & 6189 !AH comm2d, status, ierr ) 6190 !AH#endif 6191 !AH ELSE 6192 !AH workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw1) & 6193 !AH = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw1) 6194 !AH ENDIF 6195 6196 DO l = iclw+1, icrw1 6197 DO n = 0, kct 6198 6199 DO i = ifl(l), ifu(l) 6200 DO k = kfl(n), kfu(n) 6201 f(k,jbc,i) = fc(n,m,l) 6250 IF ( var == 'u' .OR. var == 'w' ) THEN 6251 ! 6252 ! Substitute the necessary parentgrid data to the work array workarrc_sn. 6253 workarrc_sn = 0.0_wp 6254 IF ( pdims(1) > 1 ) THEN 6255 #if defined( __parallel ) 6256 IF ( nxl == 0 ) THEN ! if ( bc_dirichlet_l ) 6257 workarrc_sn(0:cg%nz+1,0:2,iclw:icrw1) & 6258 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw:icrw1) 6259 ELSE IF ( nxr == nx ) THEN ! if ( bc_dirichlet_r ) 6260 workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw) & 6261 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw) 6262 ELSE 6263 workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw1) & 6264 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw1) 6265 ENDIF 6266 ! 6267 ! Leftright exchange if more than one PE subdomain in the xdirection. 6268 ! Note that in case of 3D nesting the left (pleft == MPI_PROC_NULL) and 6269 ! right (pright == MPI_PROC_NULL) boundaries are not exchanged because 6270 ! the nest domain is not cyclic. 6271 ! From left to right 6272 CALL MPI_SENDRECV( workarrc_sn(0,0,iclw+1), 1, & 6273 workarrc_sn_exchange_type, pleft, 0, & 6274 workarrc_sn(0,0,icrw), 1, & 6275 workarrc_sn_exchange_type, pright, 0, & 6276 comm2d, status, ierr ) 6277 ! 6278 ! From right to left 6279 CALL MPI_SENDRECV( workarrc_sn(0,0,icrw1), 1, & 6280 workarrc_sn_exchange_type, pright, 1, & 6281 workarrc_sn(0,0,iclw), 1, & 6282 workarrc_sn_exchange_type, pleft, 1, & 6283 comm2d, status, ierr ) 6284 #endif 6285 ELSE 6286 workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw1) & 6287 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw1) 6288 ENDIF 6289 6290 IF ( var == 'u' ) THEN 6291 6292 DO l = iclw+1, icrw 6293 DO n = 0, kct 6294 ! 6295 ! First substitute only the matchingnode values 6296 f(kfl(n):kfu(n),jbc,ifl(l)) = workarrc_sn(n,mw,l) 6297 6298 ENDDO 6299 ENDDO 6300 6301 DO l = iclw+1, icrw1 6302 DO n = 0, kct 6303 ! 6304 ! Then fill up the nodes in between with the averages 6305 DO i = ifl(l)+1, ifl(l+1)1 6306 f(kfl(n):kfu(n),jbc,i) = 0.5_wp * ( f(kfl(n):kfu(n),jbc,ifl(l)) & 6307 + f(kfl(n):kfu(n),jbc,ifl(l+1)) ) 6202 6308 ENDDO 6309 6203 6310 ENDDO 6204 6205 6311 ENDDO 6206 ENDDO 6207 6208 ELSE IF ( var == 'u' ) THEN 6209 6210 DO l = iclw+1, icrw1 6211 DO n = 0, kct 6212 6213 DO i = ifl(l), ifl(l+1)1 6214 DO k = kfl(n), kfu(n) 6215 f(k,jbc,i) = fc(n,m,l) 6216 !AH 6217 ! write(9,"('sn u: ',6(i3,2x),2(e12.5,2x))") n, m, l, k, jbc, i, fc(n,m,l), f(k,jbc,i) 6218 ! flush(9) 6219 !AH 6312 6313 ELSE IF ( var == 'w' ) THEN 6314 6315 DO l = iclw+1, icrw1 6316 DO n = 0, kct + 1 ! It is important to go up to kct+1 6317 ! 6318 ! First substitute only the matchingnode values 6319 f(kfu(n),jbc,ifl(l):ifu(l)) = workarrc_sn(n,mw,l) 6320 6321 ENDDO 6322 ENDDO 6323 6324 DO l = iclw+1, icrw1 6325 DO n = 1, kct + 1 ! It is important to go up to kct+1 6326 ! 6327 ! Then fill up the nodes in between with the averages 6328 DO k = kfu(n1)+1, kfu(n)1 6329 f(k,jbc,ifl(l):ifu(l)) = 0.5_wp * ( f(kfu(n1),jbc,ifl(l):ifu(l)) & 6330 + f(kfu(n),jbc,ifl(l):ifu(l)) ) 6220 6331 ENDDO 6332 6221 6333 ENDDO 6222 6223 6334 ENDDO 6224 ENDDO 6225 6226 ELSE IF ( var == 'w' ) THEN 6227 6228 DO l = iclw+1, icrw1 6229 DO n = 1, kct + 1 ! It is important to go up to kct+1 6230 6231 DO i = ifl(l), ifu(l) 6232 f(nzb,jbc,i) = 0.0_wp ! Because the nloop starts from n=1 instead of 0 6233 DO k = kfu(n1)+1, kfu(n) 6234 f(k,jbc,i) = fc(n,m,l) 6235 ENDDO 6236 ENDDO 6237 6238 ENDDO 6239 ENDDO 6240 6241 ELSE ! scalars 6335 6336 ENDIF 6337 6338 ELSE ! v or scalars 6242 6339 6243 6340 DO l = iclw+1, icrw1 … … 6358 6455 REAL(wp) :: fkpjp !< 6359 6456 REAL(wp) :: fkp !< 6457 REAL(wp) :: f_interp !< 6360 6458 REAL(wp) :: rcorr !< 6361 6459 REAL(wp) :: rcorr_ijk !< 6460 6461 real(wp), parameter :: c1 = 2.0_wp / 6.0_wp 6462 real(wp), parameter :: c2 = 5.0_wp / 6.0_wp 6463 real(wp), parameter :: c3 = 1.0_wp / 6.0_wp 6362 6464 6363 6465 … … 6380 6482 ENDIF 6381 6483 n = kc(k) + noff 6382 nw = 1 6383 6384 IF ( var == 'w' ) THEN 6385 !AH! 6386 !AH! Substitute the necessary parentgrid data to the work array. 6387 !AH! Note that the dimension of workarrc_t is (0:2,jcsw:jcnw,iclw:icrw), 6388 !AH! And the jc?w and ic?windex bounds depend on the location of the PE 6389 !AH! subdomain relative to the side boundaries. 6390 !AH iclc = iclw + 1 6391 !AH icrc = icrw  1 6392 !AH jcsc = jcsw + 1 6393 !AH jcnc = jcnw  1 6394 !AH IF ( bc_dirichlet_l ) THEN 6395 !AH iclc = iclw 6396 !AH ENDIF 6397 !AH IF ( bc_dirichlet_r ) THEN 6398 !AH icrc = icrw 6399 !AH ENDIF 6400 !AH IF ( bc_dirichlet_s ) THEN 6401 !AH jcsc = jcsw 6402 !AH ENDIF 6403 !AH IF ( bc_dirichlet_n ) THEN 6404 !AH jcnc = jcnw 6405 !AH ENDIF 6406 !AH workarrc_t = 0.0_wp 6407 !AH workarrc_t(0:2,jcsc:jcnc,iclc:icrc) & 6408 !AH = fc(kc(k):kc(k)+2,jcsc:jcnc,iclc:icrc) 6409 !AH! 6410 !AH! Leftright exchange if more than one PE subdomain in the xdirection. 6411 !AH! Note that in case of 3D nesting the left and right boundaries are 6412 !AH! not exchanged because the nest domain is not cyclic. 6413 !AH#if defined( __parallel ) 6414 !AH IF ( pdims(1) > 1 ) THEN 6415 !AH! 6416 !AH! From left to right 6417 !AH CALL MPI_SENDRECV( workarrc_t(0,jcsw,iclw+1), 1, & 6418 !AH workarrc_t_exchange_type_y, pleft, 0, & 6419 !AH workarrc_t(0,jcsw,icrw), 1, & 6420 !AH workarrc_t_exchange_type_y, pright, 0, & 6421 !AH comm2d, status, ierr ) 6422 !AH! 6423 !AH! From right to left 6424 !AH CALL MPI_SENDRECV( workarrc_t(0,jcsw,icrw1), 1, & 6425 !AH workarrc_t_exchange_type_y, pright, 1, & 6426 !AH workarrc_t(0,jcsw,iclw), 1, & 6427 !AH workarrc_t_exchange_type_y, pleft, 1, & 6428 !AH comm2d, status, ierr ) 6429 !AH ENDIF 6430 !AH! 6431 !AH! Southnorth exchange if more than one PE subdomain in the ydirection. 6432 !AH! Note that in case of 3D nesting the south and north boundaries are 6433 !AH! not exchanged because the nest domain is not cyclic. 6434 !AH IF ( pdims(2) > 1 ) THEN 6435 !AH! 6436 !AH! From south to north 6437 !AH CALL MPI_SENDRECV( workarrc_t(0,jcsw+1,iclw), 1, & 6438 !AH workarrc_t_exchange_type_x, psouth, 2, & 6439 !AH workarrc_t(0,jcnw,iclw), 1, & 6440 !AH workarrc_t_exchange_type_x, pnorth, 2, & 6441 !AH comm2d, status, ierr ) 6442 !AH! 6443 !AH! From north to south 6444 !AH CALL MPI_SENDRECV( workarrc_t(0,jcnw1,iclw), 1, & 6445 !AH workarrc_t_exchange_type_x, pnorth, 3, & 6446 !AH workarrc_t(0,jcsw,iclw), 1, & 6447 !AH workarrc_t_exchange_type_x, psouth, 3, & 6448 !AH comm2d, status, ierr ) 6449 !AH ENDIF 6450 !AH#endif 6451 !AH 6452 6484 nw = noff 6485 ! write(9,"(a1,2x,5(i3,2x))") var, k, kc(k), noff, n, nw 6486 ! flush(9) 6487 ! 6488 ! Substitute the necessary parentgrid data to the work array. 6489 ! Note that the dimension of workarrc_t is (0:2,jcsw:jcnw,iclw:icrw), 6490 ! And the jc?w and ic?windex bounds depend on the location of the PE 6491 ! subdomain relative to the side boundaries. 6492 iclc = iclw + 1 6493 icrc = icrw  1 6494 jcsc = jcsw + 1 6495 jcnc = jcnw  1 6496 IF ( bc_dirichlet_l ) THEN 6497 iclc = iclw 6498 ENDIF 6499 IF ( bc_dirichlet_r ) THEN 6500 icrc = icrw 6501 ENDIF 6502 IF ( bc_dirichlet_s ) THEN 6503 jcsc = jcsw 6504 ENDIF 6505 IF ( bc_dirichlet_n ) THEN 6506 jcnc = jcnw 6507 ENDIF 6508 workarrc_t = 0.0_wp 6509 workarrc_t(0:2,jcsc:jcnc,iclc:icrc) = fc(kc(k):kc(k)+2,jcsc:jcnc,iclc:icrc) 6510 ! 6511 ! Leftright exchange if more than one PE subdomain in the xdirection. 6512 ! Note that in case of 3D nesting the left and right boundaries are 6513 ! not exchanged because the nest domain is not cyclic. 6514 #if defined( __parallel ) 6515 IF ( pdims(1) > 1 ) THEN 6516 ! 6517 ! From left to right 6518 CALL MPI_SENDRECV( workarrc_t(0,jcsw,iclw+1), 1, & 6519 workarrc_t_exchange_type_y, pleft, 0, & 6520 workarrc_t(0,jcsw,icrw), 1, & 6521 workarrc_t_exchange_type_y, pright, 0, & 6522 comm2d, status, ierr ) 6523 ! 6524 ! From right to left 6525 CALL MPI_SENDRECV( workarrc_t(0,jcsw,icrw1), 1, & 6526 workarrc_t_exchange_type_y, pright, 1, & 6527 workarrc_t(0,jcsw,iclw), 1, & 6528 workarrc_t_exchange_type_y, pleft, 1, & 6529 comm2d, status, ierr ) 6530 ENDIF 6531 ! 6532 ! Southnorth exchange if more than one PE subdomain in the ydirection. 6533 ! Note that in case of 3D nesting the south and north boundaries are 6534 ! not exchanged because the nest domain is not cyclic. 6535 IF ( pdims(2) > 1 ) THEN 6536 ! 6537 ! From south to north 6538 CALL MPI_SENDRECV( workarrc_t(0,jcsw+1,iclw), 1, & 6539 workarrc_t_exchange_type_x, psouth, 2, & 6540 workarrc_t(0,jcnw,iclw), 1, & 6541 workarrc_t_exchange_type_x, pnorth, 2, & 6542 comm2d, status, ierr ) 6543 ! 6544 ! From north to south 6545 CALL MPI_SENDRECV( workarrc_t(0,jcnw1,iclw), 1, & 6546 workarrc_t_exchange_type_x, pnorth, 3, & 6547 workarrc_t(0,jcsw,iclw), 1, & 6548 workarrc_t_exchange_type_x, psouth, 3, & 6549 comm2d, status, ierr ) 6550 ENDIF 6551 #endif 6552 6553 IF ( var == 'w' ) THEN 6554 ! DO l = iclw+1, icrw1 6555 ! DO m = jcsw+1, jcnw1 6556 DO l = iclw, icrw ! Also fill up beyond edges 6557 DO m = jcsw, jcnw ! Also fill up beyond edges 6558 6559 DO i = ifl(l), ifu(l) 6560 DO j = jfl(m), jfu(m) 6561 f(k,j,i) = workarrc_t(nw,m,l) 6562 ENDDO 6563 ENDDO 6564 6565 ENDDO 6566 ENDDO 6567 ENDIF 6568 6569 IF ( var == 'u' ) THEN 6570 6571 DO l = iclw+1, icrw ! Note that we have to go up to icrw here. 6572 DO m = jcsw+1, jcnw1 6573 ! 6574 ! 1storder upwind interpolation 6575 f_interp = workarrc_t(nw,m,l) 6576 ! 6577 ! 3rdorder upwind biased interpolation. Note that w is available 6578 ! here since it is already interpolated. 6579 ! IF ( w(k,jfl(m),ifl(l)) + w(k,jfl(m),ifl(l)1) < 0.0_wp ) THEN 6580 ! f_interp = c1 * workarrc_t(0,m,l) + c2 * workarrc_t(1,m,l) & 6581 ! + c3 * workarrc_t(2,m,l) 6582 ! ELSE 6583 ! f_interp = workarrc_t(nw,m,l) 6584 ! ENDIF 6585 6586 ! DO i = ifl(l), ifl(l+1)1 6587 ! DO j = jfl(m), jfu(m) 6588 ! f(k,j,i) = fc(n,m,l) 6589 ! ENDDO 6590 ! ENDDO 6591 ! 6592 ! First substitute only the matchingnode values 6593 f(k,jfl(m):jfu(m),ifl(l)) = f_interp 6594 6595 ENDDO 6596 ENDDO 6597 ! 6598 ! Then fill up the nodes in between with the averages 6453 6599 DO l = iclw+1, icrw1 6454 6600 DO m = jcsw+1, jcnw1 6455 6601 6602 DO i = ifl(l)+1, ifl(l+1)1 6603 f(k,jfl(m):jfu(m),i) = 0.5_wp * ( f(k,jfl(m):jfu(m),ifl(l)) & 6604 + f(k,jfl(m):jfu(m),ifl(l+1)) ) 6605 ENDDO 6606 6607 ENDDO 6608 ENDDO 6609 6610 ENDIF 6611 6612 IF ( var == 'v' ) THEN 6613 6614 DO l = iclw+1, icrw1 6615 DO m = jcsw+1, jcnw ! Note that we have to go up to jcnw here. 6616 ! 6617 ! 1storder upwind interpolation 6618 f_interp = workarrc_t(nw,m,l) 6619 ! 6620 ! 3rdorder upwind biased interpolation. Note that w is available 6621 ! here since it is already interpolated. 6622 ! IF ( w(k,jfl(m),ifl(l)) + w(k,jfl(m)1,ifl(l)) < 0.0_wp ) THEN 6623 ! f_interp = c1 * workarrc_t(0,m,l) + c2 * workarrc_t(1,m,l) + c3 * workarrc_t(2,m,l) 6624 ! ELSE 6625 ! f_interp = workarrc_t(nw,m,l) 6626 ! ENDIF 6627 6628 ! DO i = ifl(l), ifu(l) 6629 ! DO j = jfl(m), jfl(m+1)1 6630 ! f(k,j,i) = fc(n,m,l) 6631 ! ENDDO 6632 ! ENDDO 6633 ! 6634 ! First substitute only the matchingnode values 6635 f(k,jfl(m),ifl(l):ifu(l)) = f_interp 6636 6637 ENDDO 6638 ENDDO 6639 ! 6640 ! Then fill up the nodes in between with the averages 6641 DO l = iclw+1, icrw1 6642 DO m = jcsw+1, jcnw1 6643 6644 DO j = jfl(m)+1, jfl(m+1)1 6645 f(k,j,ifl(l):ifu(l)) = 0.5_wp * ( f(k,jfl(m),ifl(l):ifu(l)) & 6646 + f(k,jfl(m+1),ifl(l):ifu(l)) ) 6647 ENDDO 6648 6649 ENDDO 6650 ENDDO 6651 6652 ENDIF 6653 6654 IF ( var == 's' ) THEN 6655 DO l = iclw+1, icrw1 6656 DO m = jcsw+1, jcnw1 6657 ! 6658 ! 1storder upwind interpolation 6659 f_interp = workarrc_t(nw,m,l) 6660 ! 6661 ! 3rdorder upwind biased interpolation. Note that w is available 6662 ! here since it is already interpolated. 6663 ! IF ( w(k,jfl(m),ifl(l)) < 0.0_wp ) THEN 6664 ! f_interp = c1 * workarrc_t(0,m,l) + c2 * workarrc_t(1,m,l) + c3 * workarrc_t(2,m,l) 6665 ! ELSE 6666 ! f_interp = workarrc_t(nw,m,l) 6667 ! ENDIF 6668 6456 6669 DO i = ifl(l), ifu(l) 6457 6670 DO j = jfl(m), jfu(m) 6458 f(k,j,i) = f c(n,m,l)6671 f(k,j,i) = f_interp 6459 6672 ENDDO 6460 6673 ENDDO … … 6462 6675 ENDDO 6463 6676 ENDDO 6464 6465 ELSE IF ( var == 'u' ) THEN 6466 6467 DO l = iclw+1, icrw1 6468 DO m = jcsw+1, jcnw1 6469 6470 DO i = ifl(l), ifl(l+1)1 6471 DO j = jfl(m), jfu(m) 6472 f(k,j,i) = fc(n,m,l) 6473 ENDDO 6474 ENDDO 6475 6476 ENDDO 6477 ENDDO 6478 6479 ELSE IF ( var == 'v' ) THEN 6480 6481 DO l = iclw+1, icrw1 6482 DO m = jcsw+1, jcnw1 6483 6484 DO i = ifl(l), ifu(l) 6485 DO j = jfl(m), jfl(m+1)1 6486 f(k,j,i) = fc(n,m,l) 6487 ENDDO 6488 ENDDO 6489 6490 ENDDO 6491 ENDDO 6492 6493 ELSE ! scalars 6494 6495 DO l = iclw+1, icrw1 6496 DO m = jcsw+1, jcnw1 6497 6498 DO i = ifl(l), ifu(l) 6499 DO j = jfl(m), jfu(m) 6500 f(k,j,i) = fc(n,m,l) 6501 ENDDO 6502 ENDDO 6503 6504 ENDDO 6505 ENDDO 6506 6507 ENDIF ! var 6677 ENDIF 6508 6678 ! 6509 6679 ! Just fill up the redundant second ghostnode layer in case of var == w.
Note: See TracChangeset
for help on using the changeset viewer.