Changeset 3741 for palm/trunk/SOURCE
- Timestamp:
- Feb 13, 2019 4:24:49 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/pmc_interface_mod.f90 (modified) (39 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/pmc_interface_mod.f90
r3708 r3741 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Interpolations and child initialization adjusted to handle set ups with child 28 ! pe-subdomain dimension not integer divisible by the grid-spacing ratio in the 29 ! respective direction. Set ups with pe-subdomain dimension smaller than the 30 ! grid-spacing ratio in the respective direction are now forbidden. 31 ! 32 ! 3708 2019-01-30 12:58:13Z hellstea 27 33 ! Checks for parent / child grid line matching introduced. 28 34 ! Interpolation of nest-boundary-tangential velocity components revised. … … 1504 1510 INTEGER(iwp) :: jauxs !< 1505 1511 INTEGER(iwp) :: jauxn !< 1506 REAL(wp) :: loffset !<1507 REAL(wp) :: noffset !<1508 REAL(wp) :: roffset !<1509 REAL(wp) :: soffset !<1510 1512 REAL(wp) :: xexl !< Parent-grid array exceedance behind the left edge of the child PE subdomain 1511 1513 REAL(wp) :: xexr !< Parent-grid array exceedance behind the right edge of the child PE subdomain … … 1580 1582 !-- Left 1581 1583 IF ( bc_dirichlet_l ) THEN 1582 loffset = MOD( coord_x(nxl), cg%dx ) 1583 xexl = 2 * cg%dx + loffset 1584 xexl = 2 * cg%dx 1584 1585 iauxl = 0 1585 1586 ELSE … … 1597 1598 !-- Right 1598 1599 IF ( bc_dirichlet_r ) THEN 1599 roffset = MOD( coord_x(nxr+1), cg%dx ) 1600 xexr = 2 * cg%dx + roffset 1600 xexr = 2 * cg%dx 1601 1601 iauxr = 0 1602 1602 ELSE … … 1607 1607 DO i = cg%nx, 0 , -1 1608 1608 IF ( cg%coord_x(i) + 0.5_wp * cg%dx <= xce ) THEN 1609 icr = MIN( cg%nx, i)1609 icr = MIN( cg%nx, MAX( icl, i ) ) 1610 1610 EXIT 1611 1611 ENDIF … … 1614 1614 !-- South 1615 1615 IF ( bc_dirichlet_s ) THEN 1616 soffset = MOD( coord_y(nys), cg%dy ) 1617 yexs = 2 * cg%dy + soffset 1616 yexs = 2 * cg%dy 1618 1617 jauxs = 0 1619 1618 ELSE … … 1631 1630 !-- North 1632 1631 IF ( bc_dirichlet_n ) THEN 1633 noffset = MOD( coord_y(nyn+1), cg%dy ) 1634 yexn = 2 * cg%dy + noffset 1632 yexn = 2 * cg%dy 1635 1633 jauxn = 0 1636 1634 ELSE … … 1641 1639 DO j = cg%ny, 0 , -1 1642 1640 IF ( cg%coord_y(j) + 0.5_wp * cg%dy <= yce ) THEN 1643 jcn = MIN( cg%ny, j)1641 jcn = MIN( cg%ny, MAX( jcs, j ) ) 1644 1642 EXIT 1645 1643 ENDIF … … 3556 3554 ! 3557 3555 !-- i-indices of u for each ii-index value 3558 !-- ii=icr is redundant for anterpolation3559 tolerance = 0.000001_wp * dx3560 3556 istart = nxlg 3561 !AH DO ii = icl, icr 3562 !AH DO ii = icl-1, icr+1 3563 DO ii = icla, icra 3564 3565 ! 3566 !-- In case the child and parent grid lines match in x 3567 !-- use only the local k,j-child-grid plane for the anterpolation, 3568 !-- i.e use 2-D anterpolation. Else, use a 3-D anterpolation. 3557 DO ii = icla, icra 3558 ! 3559 !-- The parent and child grid lines do always match in x, hence we 3560 !-- use only the local k,j-child-grid plane for the anterpolation. 3569 3561 i = istart 3570 DO WHILE ( coord_x(i) < cg%coord_x(ii) - tolerance .AND.i < nxrg )3562 DO WHILE ( coord_x(i) < cg%coord_x(ii) .AND. i < nxrg ) 3571 3563 i = i + 1 3572 3564 ENDDO 3573 IF ( ABS( coord_x(i) - cg%coord_x(ii) ) < tolerance ) THEN 3574 i = istart 3575 DO WHILE ( coord_x(i) < cg%coord_x(ii) .AND. i < nxrg ) 3576 i = i + 1 3577 ENDDO 3578 iflu(ii) = MIN( MAX( i, nxlg ), nxrg ) 3579 ifuu(ii) = iflu(ii) 3580 istart = iflu(ii) 3581 ELSE 3582 i = istart 3583 DO WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx ) & 3584 .AND. ( i < nxrg ) ) 3585 i = i + 1 3586 ENDDO 3587 iflu(ii) = MIN( MAX( i, nxlg ), nxrg ) 3588 ir = i 3589 DO WHILE ( ( coord_x(ir) <= cg%coord_x(ii) + 0.5_wp * cg%dx ) & 3590 .AND. ( i < nxrg+1 ) ) 3591 i = i + 1 3592 ir = MIN( i, nxrg ) 3593 ENDDO 3594 ifuu(ii) = MIN( MAX( i-1, iflu(ii) ), nxrg ) 3595 istart = iflu(ii) 3596 ENDIF 3565 iflu(ii) = MIN( MAX( i, nxlg ), nxrg ) 3566 ifuu(ii) = iflu(ii) 3567 istart = iflu(ii) 3597 3568 ! 3598 3569 !-- Print out the index bounds for checking and debugging purposes … … 3605 3576 ! 3606 3577 !-- i-indices of others for each ii-index value 3607 !-- ii=icr is redundant for anterpolation3608 3578 istart = nxlg 3609 !AH DO ii = icl, icr3610 !AH DO ii = icl-1, icr+13611 3579 DO ii = icla, icra 3612 3580 i = istart … … 3633 3601 ! 3634 3602 !-- j-indices of v for each jj-index value 3635 !-- jj=jcn is redundant for anterpolation3636 tolerance = 0.000001_wp * dy3637 3603 jstart = nysg 3638 !AH DO jj = jcs, jcn3639 !AH DO jj = jcs-1, jcn+13640 3604 DO jj = jcsa, jcna 3641 3605 ! 3642 !-- In case the child and parent grid lines match in y 3643 !-- use only the local k,i-child-grid plane for the anterpolation, 3644 !-- i.e use 2-D anterpolation. Else, use a 3-D anterpolation. 3606 !-- The parent and child grid lines do always match in y, hence we 3607 !-- use only the local k,i-child-grid plane for the anterpolation. 3645 3608 j = jstart 3646 DO WHILE ( coord_y(j) < cg%coord_y(jj) - tolerance .AND.j < nyng )3609 DO WHILE ( coord_y(j) < cg%coord_y(jj) .AND. j < nyng ) 3647 3610 j = j + 1 3648 3611 ENDDO 3649 IF ( ABS( coord_y(j) - cg%coord_y(jj) ) < tolerance ) THEN 3650 j = jstart 3651 DO WHILE ( coord_y(j) < cg%coord_y(jj) .AND. j < nyng ) 3652 j = j + 1 3653 ENDDO 3654 jflv(jj) = MIN( MAX( j, nysg ), nyng ) 3655 jfuv(jj) = jflv(jj) 3656 jstart = jflv(jj) 3657 ELSE 3658 j = jstart 3659 DO WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy ) & 3660 .AND. ( j < nyng ) ) 3661 j = j + 1 3662 ENDDO 3663 jflv(jj) = MIN( MAX( j, nysg ), nyng ) 3664 jr = j 3665 DO WHILE ( ( coord_y(jr) <= cg%coord_y(jj) + 0.5_wp * cg%dy ) & 3666 .AND. ( j < nyng+1 ) ) 3667 j = j + 1 3668 jr = MIN( j, nyng ) 3669 ENDDO 3670 jfuv(jj) = MIN( MAX( j-1, jflv(jj) ), nyng ) 3671 jstart = jflv(jj) 3672 ENDIF 3612 jflv(jj) = MIN( MAX( j, nysg ), nyng ) 3613 jfuv(jj) = jflv(jj) 3614 jstart = jflv(jj) 3673 3615 ! 3674 3616 !-- Print out the index bounds for checking and debugging purposes … … 3680 3622 ! 3681 3623 !-- j-indices of others for each jj-index value 3682 !-- jj=jcn is redundant for anterpolation3683 3624 jstart = nysg 3684 !AH DO jj = jcs, jcn3685 !AH DO jj = jcs-1, jcn+13686 3625 DO jj = jcsa, jcna 3687 3626 j = jstart … … 3709 3648 !-- k-indices of w for each kk-index value 3710 3649 !-- Note that anterpolation index limits are needed also for the top boundary 3711 !-- ghost cell level because of the reversibility correctionin the interpolation.3650 !-- ghost cell level because they are used also in the interpolation. 3712 3651 kstart = 0 3713 3652 kflw(0) = 0 3714 3653 kfuw(0) = 0 3715 tolerance = 0.000001_wp * dzw(1)3716 3654 DO kk = 1, cg%nz+1 3717 3655 ! 3718 !-- In case the child and parent grid lines match in z 3719 !-- use only the local j,i-child-grid plane for the anterpolation, 3720 !-- i.e use 2-D anterpolation. Else, use a 3-D anterpolation. 3656 !-- The parent and child grid lines do always match in z, hence we 3657 !-- use only the local j,i-child-grid plane for the anterpolation. 3721 3658 k = kstart 3722 DO WHILE ( zw(k) < cg%zw(kk) - tolerance .AND. k < nzt+1)3659 DO WHILE ( ( zw(k) < cg%zw(kk) ) .AND. ( k < nzt+1 ) ) 3723 3660 k = k + 1 3724 3661 ENDDO 3725 IF ( ABS( zw(k) - cg%zw(kk) ) < tolerance ) THEN 3726 k = kstart 3727 DO WHILE ( ( zw(k) < cg%zw(kk) ) .AND. ( k < nzt+1 ) ) 3728 k = k + 1 3729 ENDDO 3730 kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 3731 kfuw(kk) = kflw(kk) 3732 kstart = kflw(kk) 3733 ELSE 3734 k = kstart 3735 DO WHILE ( ( zw(k) < cg%zu(kk) ) .AND. ( k < nzt ) ) 3736 k = k + 1 3737 ENDDO 3738 kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 3739 IF ( kk+1 <= cg%nz+1 ) THEN 3740 DO WHILE ( ( zw(k) <= cg%zu(kk+1) ) .AND. ( k < nzt+1 ) ) 3741 k = k + 1 3742 IF ( k > nzt + 1 ) EXIT ! This EXIT is to prevent zu(k) from flowing over. 3743 ENDDO 3744 kfuw(kk) = MIN( MAX( k-1, kflw(kk) ), nzt + 1 ) 3745 ELSE 3746 kfuw(kk) = kflw(kk) 3747 ENDIF 3748 kstart = kflw(kk) 3749 ENDIF 3750 !AH 3662 kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 ) 3663 kfuw(kk) = kflw(kk) 3664 kstart = kflw(kk) 3665 ! 3666 !-- Print out the index bounds for checking and debugging purposes 3751 3667 write(9,"('pmci_init_anterp_tophat, kk, kflw, kfuw: ', 4(i4,2x), 2(e12.5,2x))") & 3752 3668 kk, kflw(kk), kfuw(kk), nzt, cg%zu(kk), cg%zw(kk) … … 3761 3677 ! 3762 3678 !-- Note that anterpolation index limits are needed also for the top boundary 3763 !-- ghost cell level because of the reversibility correctionin the interpolation.3679 !-- ghost cell level because they are used also in the interpolation. 3764 3680 DO kk = 1, cg%nz+1 3765 3681 k = kstart … … 3782 3698 kflo(cg%nz+1) = nzt+kgsr ! Obsolete 3783 3699 kfuo(cg%nz+1) = nzt+kgsr ! Obsolete 3784 3700 ! 3701 !-- Print out the index bounds for checking and debugging purposes 3785 3702 DO kk = 1, cg%nz+1 3786 3703 write(9,"('pmci_init_anterp_tophat, kk, kflo, kfuo: ', 4(i4,2x), 2(e12.5,2x))") & … … 3789 3706 ENDDO 3790 3707 write(9,*) 3791 !AH3792 3793 3708 ! 3794 3709 !-- Precomputation of number of fine-grid nodes inside parent-grid cells. … … 3910 3825 SUBROUTINE pmci_check_grid_matching 3911 3826 ! 3912 !-- Check that the grid lines of child and parent do match 3827 !-- Check that the grid lines of child and parent do match. 3828 !-- Also check that the child subdomain width is not smaller than 3829 !-- the parent grid spacing in the respective direction. 3913 3830 IMPLICIT NONE 3914 REAL(wp), PARAMETER :: tolefac = 1.0E-6_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 3831 REAL(wp), PARAMETER :: tolefac = 1.0E-6_wp !< Relative tolerence for grid-line matching 3832 REAL(wp) :: pesdwx !< Child subdomain width in x-direction 3833 REAL(wp) :: pesdwy !< Child subdomain width in y-direction 3834 REAL(wp) :: tolex !< Tolerance for grid-line matching in x-direction 3835 REAL(wp) :: toley !< Tolerance for grid-line matching in y-direction 3836 REAL(wp) :: tolez !< Tolerance for grid-line matching in z-direction 3837 INTEGER(iwp) :: non_matching_lower_left_corner !< Flag for non-matching lower left corner 3838 INTEGER(iwp) :: non_int_gsr_x !< Flag for non-integer grid-spacing ration in x-direction 3839 INTEGER(iwp) :: non_int_gsr_y !< Flag for non-integer grid-spacing ration in y-direction 3840 INTEGER(iwp) :: non_int_gsr_z !< Flag for non-integer grid-spacing ration in z-direction 3841 INTEGER(iwp) :: too_narrow_pesd_x !< Flag for too narrow pe-subdomain in x-direction 3842 INTEGER(iwp) :: too_narrow_pesd_y !< Flag for too narrow pe-subdomain in y-direction 3922 3843 3923 3844 … … 3926 3847 non_int_gsr_y = 0 3927 3848 non_int_gsr_z = 0 3849 too_narrow_pesd_x = 0 3850 too_narrow_pesd_y = 0 3928 3851 3929 3852 IF ( myid == 0 ) THEN … … 3944 3867 ENDDO 3945 3868 3869 pesdwx = REAL( nxr - nxl + 1, KIND=wp ) 3870 IF ( pesdwx / REAL( igsr, KIND=wp ) < 1.0_wp ) too_narrow_pesd_x = 1 3871 pesdwy = REAL( nyn - nys + 1, KIND=wp ) 3872 IF ( pesdwy / REAL( jgsr, KIND=wp ) < 1.0_wp ) too_narrow_pesd_y = 1 3873 3946 3874 ENDIF 3947 3875 … … 3973 3901 '( parent dz / child dz ) must have an integer value for each z-level' 3974 3902 CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 ) 3975 ENDIF 3903 ENDIF 3904 3905 CALL MPI_BCAST( too_narrow_pesd_x, 1, MPI_INTEGER, 0, comm2d, ierr ) 3906 IF ( too_narrow_pesd_x > 0 ) THEN 3907 WRITE ( message_string, * ) 'child subdomain width in x-direction ', & 3908 'must not be smaller than its parent grid dx. Change the PE-grid ', & 3909 'setting (npex, npey) to satisfy this requirement.' 3910 CALL message( 'pmci_check_grid_matching', 'PA0587', 3, 2, 0, 6, 0 ) 3911 ENDIF 3912 3913 CALL MPI_BCAST( too_narrow_pesd_y, 1, MPI_INTEGER, 0, comm2d, ierr ) 3914 IF ( too_narrow_pesd_y > 0 ) THEN 3915 WRITE ( message_string, * ) 'child subdomain width in y-direction ', & 3916 'must not be smaller than its parent grid dy. Change the PE-grid ', & 3917 'setting (npex, npey) to satisfy this requirement.' 3918 CALL message( 'pmci_check_grid_matching', 'PA0587', 3, 2, 0, 6, 0 ) 3919 ENDIF 3976 3920 3977 3921 END SUBROUTINE pmci_check_grid_matching … … 4521 4465 ! 4522 4466 !-- Interpolation of the internal values for the child-domain initialization 4523 !-- This subroutine is based on trilinear interpolation.4524 !-- Coding based on interp_tril_lr/sn/t4525 4467 IMPLICIT NONE 4526 4468 … … 4558 4500 4559 4501 INTEGER(iwp) :: i !< 4502 INTEGER(iwp) :: ib !< 4503 INTEGER(iwp) :: ie !< 4504 INTEGER(iwp) :: ifirst !< 4505 INTEGER(iwp) :: ilast !< 4560 4506 INTEGER(iwp) :: j !< 4507 INTEGER(iwp) :: jb !< 4508 INTEGER(iwp) :: je !< 4509 INTEGER(iwp) :: jfirst !< 4510 INTEGER(iwp) :: jlast !< 4561 4511 INTEGER(iwp) :: k !< 4562 4512 INTEGER(iwp) :: l !< … … 4574 4524 mb = jcs 4575 4525 me = jcn 4526 ifirst = nxl 4527 ilast = nxr 4528 jfirst = nys 4529 jlast = nyn 4576 4530 4577 4531 IF ( nesting_mode /= 'vertical' ) THEN 4578 4532 IF ( bc_dirichlet_l ) THEN 4579 lb = icl + 1 4533 lb = icl + 1 4534 ifirst = nxl - 1 4580 4535 ! 4581 4536 !-- For u, nxl is a ghost node, but not for the other variables 4582 4537 IF ( var == 'u' ) THEN 4583 lb = icl + 2 4538 lb = icl + 2 4539 ifirst = nxl 4584 4540 ENDIF 4585 4541 ENDIF 4586 4542 IF ( bc_dirichlet_s ) THEN 4587 mb = jcs + 1 4543 mb = jcs + 1 4544 jfirst = nys - 1 4588 4545 ! 4589 4546 !-- For v, nys is a ghost node, but not for the other variables 4590 4547 IF ( var == 'v' ) THEN 4591 mb = jcs + 2 4548 mb = jcs + 2 4549 jfirst = nys 4592 4550 ENDIF 4593 4551 ENDIF 4594 4552 IF ( bc_dirichlet_r ) THEN 4595 le = icr - 1 4553 le = icr - 1 4554 ilast = nxr + 1 4596 4555 ENDIF 4597 4556 IF ( bc_dirichlet_n ) THEN 4598 me = jcn - 1 4557 me = jcn - 1 4558 jlast = nyn + 1 4599 4559 ENDIF 4600 4560 ENDIF … … 4615 4575 IF ( var == 'u' ) THEN 4616 4576 4577 ib = ifl(lb) 4578 ie = ifl(le+1) - 1 4579 jb = jfl(mb) 4580 je = jfu(me) 4617 4581 DO l = lb, le 4618 4582 DO m = mb, me … … 4633 4597 ELSE IF ( var == 'v' ) THEN 4634 4598 4599 ib = ifl(lb) 4600 ie = ifu(le) 4601 jb = jfl(mb) 4602 je = jfl(me+1) - 1 4635 4603 DO l = lb, le 4636 4604 DO m = mb, me … … 4651 4619 ELSE IF ( var == 'w' ) THEN 4652 4620 4621 ib = ifl(lb) 4622 ie = ifu(le) 4623 jb = jfl(mb) 4624 je = jfu(me) 4653 4625 DO l = lb, le 4654 4626 DO m = mb, me … … 4670 4642 ELSE ! scalars 4671 4643 4644 ib = ifl(lb) 4645 ie = ifu(le) 4646 jb = jfl(mb) 4647 je = jfu(me) 4672 4648 DO l = lb, le 4673 4649 DO m = mb, me … … 4687 4663 4688 4664 ENDIF ! var 4665 ! 4666 !-- If the subdomain i- and/or j-dimension (nx/npex and/or ny/npey) is 4667 !-- not integer divisible by the grid spacing ratio in its direction, 4668 !-- the above loops will return with unfilled gaps in the initial fields. 4669 !-- These gaps, if present, are filled here. 4670 IF ( ib > ifirst ) THEN 4671 DO i = ifirst, ib-1 4672 f(:,:,i) = f(:,:,ib) 4673 ENDDO 4674 ENDIF 4675 IF ( ie < ilast ) THEN 4676 DO i = ie+1, ilast 4677 f(:,:,i) = f(:,:,ie) 4678 ENDDO 4679 ENDIF 4680 IF ( jb > jfirst ) THEN 4681 DO j = jfirst, jb-1 4682 f(:,j,:) = f(:,jb,:) 4683 ENDDO 4684 ENDIF 4685 IF ( je < jlast ) THEN 4686 DO j = je+1, jlast 4687 f(:,j,:) = f(:,je,:) 4688 ENDDO 4689 ENDIF 4689 4690 4690 4691 END SUBROUTINE pmci_interp_1sto_all … … 5754 5755 ! 5755 5756 !-- Interpolation of ghost-node values used as the child-domain boundary 5756 !-- conditions. This subroutine handles the left and right boundaries. It is 5757 !-- based on trilinear interpolation. 5758 5757 !-- conditions. This subroutine handles the left and right boundaries. 5759 5758 IMPLICIT NONE 5760 5759 … … 5865 5864 REAL(wp) :: rcorr_ijk !< Reversibility correction distributed to the individual child-grid nodes 5866 5865 5867 real(wp), parameter :: c1 = 2.0_wp / 6.0_wp5868 real(wp), parameter :: c2 = 5.0_wp / 6.0_wp5869 real(wp), parameter :: c3 = -1.0_wp / 6.0_wp5866 ! real(wp), parameter :: c1 = 2.0_wp / 6.0_wp 5867 ! real(wp), parameter :: c2 = 5.0_wp / 6.0_wp 5868 ! real(wp), parameter :: c3 = -1.0_wp / 6.0_wp 5870 5869 5871 5870 ! … … 5940 5939 var_flag = 0 5941 5940 ENDIF 5942 5943 IF ( var == 'v' .OR. var == 'w' ) THEN 5944 ! 5945 !-- Substitute the necessary parent-grid data to the work array workarrc_lr. 5946 workarrc_lr = 0.0_wp 5947 IF ( pdims(2) > 1 ) THEN 5941 ! 5942 !-- Substitute the necessary parent-grid data to the work array workarrc_lr. 5943 workarrc_lr = 0.0_wp 5944 IF ( pdims(2) > 1 ) THEN 5948 5945 #if defined( __parallel ) 5949 IF ( nys == 0 ) THEN 5950 workarrc_lr(0:cg%nz+1,jcsw:jcnw-1,0:2) & 5951 = fc(0:cg%nz+1,jcsw:jcnw-1,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:jcnw-1,0:2) & 5957 = fc(0:cg%nz+1,jcsw+1:jcnw-1,lbeg:lbeg+2) 5958 ENDIF 5959 ! 5960 !-- South-north exchange if more than one PE subdomain in the y-direction. 5961 !-- Note that in case of 3-D 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,jcnw-1,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 5946 IF ( bc_dirichlet_s ) THEN 5947 workarrc_lr(0:cg%nz+1,jcsw:jcnw-1,0:2) & 5948 = fc(0:cg%nz+1,jcsw:jcnw-1,lbeg:lbeg+2) 5949 ELSE IF ( bc_dirichlet_n ) THEN 5950 workarrc_lr(0:cg%nz+1,jcsw+1:jcnw,0:2) & 5951 = fc(0:cg%nz+1,jcsw+1:jcnw,lbeg:lbeg+2) 5978 5952 ELSE 5979 workarrc_lr(0:cg%nz+1,jcsw :jcnw,0:2)&5980 = fc(0:cg%nz+1,jcsw :jcnw,lbeg:lbeg+2)5953 workarrc_lr(0:cg%nz+1,jcsw+1:jcnw-1,0:2) & 5954 = fc(0:cg%nz+1,jcsw+1:jcnw-1,lbeg:lbeg+2) 5981 5955 ENDIF 5982 5983 IF ( var == 'v' ) THEN 5956 ! 5957 !-- South-north exchange if more than one PE subdomain in the y-direction. 5958 !-- Note that in case of 3-D nesting the south (psouth == MPI_PROC_NULL) 5959 !-- and north (pnorth == MPI_PROC_NULL) boundaries are not exchanged 5960 !-- because the nest domain is not cyclic. 5961 !-- From south to north 5962 CALL MPI_SENDRECV( workarrc_lr(0,jcsw+1,0), 1, & 5963 workarrc_lr_exchange_type, psouth, 0, & 5964 workarrc_lr(0,jcnw,0), 1, & 5965 workarrc_lr_exchange_type, pnorth, 0, & 5966 comm2d, status, ierr ) 5967 ! 5968 !-- From north to south 5969 CALL MPI_SENDRECV( workarrc_lr(0,jcnw-1,0), 1, & 5970 workarrc_lr_exchange_type, pnorth, 1, & 5971 workarrc_lr(0,jcsw,0), 1, & 5972 workarrc_lr_exchange_type, psouth, 1, & 5973 comm2d, status, ierr ) 5974 #endif 5975 ELSE 5976 workarrc_lr(0:cg%nz+1,jcsw:jcnw,0:2) & 5977 = fc(0:cg%nz+1,jcsw:jcnw,lbeg:lbeg+2) 5978 ENDIF 5979 5980 IF ( var == 'v' ) THEN 5984 5981 5985 DO m = jcsw+1, jcnw 5986 DO n = 0, kct 5987 ! 5988 !-- First substitute only the matching-node values 5989 f(kfl(n):kfu(n),jfl(m),ibc) = workarrc_lr(n,m,lw) 5990 5982 DO m = jcsw, jcnw-1 5983 DO n = 0, kct 5984 ! 5985 !-- Use averages of the neighbouring matching grid-line values 5986 DO j = jfl(m), jfl(m+1) 5987 f(kfl(n):kfu(n),j,ibc) = 0.5_wp * ( workarrc_lr(n,m,lw) & 5988 + workarrc_lr(n,m+1,lw) ) 5989 ENDDO 5990 ! 5991 !-- Then set the values along the matching grid-lines 5992 IF ( MOD( jfl(m), jgsr ) == 0 ) THEN 5993 f(kfl(n):kfu(n),jfl(m),ibc) = workarrc_lr(n,m,lw) 5994 ENDIF 5995 ENDDO 5996 ENDDO 5997 ! 5998 !-- Finally, set the values along the last matching grid-line 5999 IF ( MOD( jfl(jcnw), jgsr ) == 0 ) THEN 6000 DO n = 0, kct 6001 f(kfl(n):kfu(n),jfl(jcnw),ibc) = workarrc_lr(n,jcnw,lw) 6002 ENDDO 6003 ENDIF 6004 ! 6005 !-- A gap may still remain in some cases if the subdomain size is not 6006 !-- divisible by the grid-spacing ratio. In such a case, fill the 6007 !-- gap. Note however, this operation may produce some additional 6008 !-- momentum conservation error. 6009 IF ( jfl(jcnw) < nyn ) THEN 6010 DO n = 0, kct 6011 DO j = jfl(jcnw)+1, nyn 6012 f(kfl(n):kfu(n),j,ibc) = f(kfl(n):kfu(n),jfl(jcnw),ibc) 5991 6013 ENDDO 5992 6014 ENDDO 5993 5994 DO m = jcsw+1, jcnw-1 5995 DO n = 0, kct5996 ! 5997 !-- Then fill up the nodes in between with the averages 5998 DO j = jfl(m)+1, jfl(m+1)-15999 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) ) 6001 ENDDO6015 ENDIF 6016 6017 ELSE IF ( var == 'w' ) THEN 6018 6019 DO m = jcsw, jcnw 6020 DO n = 0, kct + 1 ! It is important to go up to kct+1 6021 ! 6022 !-- First substitute only the matching-node values 6023 f(kfu(n),jfl(m):jfu(m),ibc) = workarrc_lr(n,m,lw) 6002 6024 6025 ENDDO 6026 ENDDO 6027 6028 DO m = jcsw, jcnw 6029 DO n = 1, kct + 1 ! It is important to go up to kct+1 6030 ! 6031 !-- Then fill up the nodes in between with the averages 6032 DO k = kfu(n-1)+1, kfu(n)-1 6033 f(k,jfl(m):jfu(m),ibc) = 0.5_wp * ( f(kfu(n-1),jfl(m):jfu(m),ibc) & 6034 + f(kfu(n),jfl(m):jfu(m),ibc) ) 6003 6035 ENDDO 6036 6004 6037 ENDDO 6005 6006 ELSE IF ( var == 'w' ) THEN 6007 6008 DO m = jcsw+1, jcnw-1 6009 DO n = 0, kct + 1 ! It is important to go up to kct+1 6010 ! 6011 !-- First substitute only the matching-node 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, jcnw-1 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(n-1)+1, kfu(n)-1 6022 f(k,jfl(m):jfu(m),ibc) = 0.5_wp * ( f(kfu(n-1),jfl(m):jfu(m),ibc) & 6023 + f(kfu(n),jfl(m):jfu(m),ibc) ) 6024 ENDDO 6025 6026 ENDDO 6027 ENDDO 6028 6029 ENDIF 6038 ENDDO 6030 6039 6031 6040 ELSE ! u or scalars 6032 6041 6033 DO m = jcsw +1, jcnw-16042 DO m = jcsw, jcnw 6034 6043 DO n = 0, kct 6035 6044 6036 6045 DO j = jfl(m), jfu(m) 6037 6046 DO k = kfl(n), kfu(n) 6038 f(k,j,ibc) = fc(n,m,l)6047 f(k,j,ibc) = workarrc_lr(n,m,lw) 6039 6048 ENDDO 6040 6049 ENDDO … … 6069 6078 !-- Interpolation of ghost-node values used as the child-domain boundary 6070 6079 !-- conditions. This subroutine handles the south and north boundaries. 6071 !-- This subroutine is based on trilinear interpolation.6072 6073 6080 IMPLICIT NONE 6074 6081 … … 6172 6179 REAL(wp) :: rcorr_ijk !< Reversibility correction distributed to the individual child-grid nodes 6173 6180 6174 real(wp), parameter :: c1 = 2.0_wp / 6.0_wp6175 real(wp), parameter :: c2 = 5.0_wp / 6.0_wp6176 real(wp), parameter :: c3 = -1.0_wp / 6.0_wp6181 ! real(wp), parameter :: c1 = 2.0_wp / 6.0_wp 6182 ! real(wp), parameter :: c2 = 5.0_wp / 6.0_wp 6183 ! real(wp), parameter :: c3 = -1.0_wp / 6.0_wp 6177 6184 6178 6185 ! … … 6247 6254 var_flag = 0 6248 6255 ENDIF 6249 6250 IF ( var == 'u' .OR. var == 'w' ) THEN 6251 ! 6252 !-- Substitute the necessary parent-grid data to the work array workarrc_sn. 6253 workarrc_sn = 0.0_wp 6254 IF ( pdims(1) > 1 ) THEN 6256 ! 6257 !-- Substitute the necessary parent-grid data to the work array workarrc_sn. 6258 workarrc_sn = 0.0_wp 6259 IF ( pdims(1) > 1 ) THEN 6255 6260 #if defined( __parallel ) 6256 IF ( nxl == 0 ) THEN ! if ( bc_dirichlet_l ) 6257 workarrc_sn(0:cg%nz+1,0:2,iclw:icrw-1) & 6258 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw:icrw-1) 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:icrw-1) & 6264 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1) 6265 ENDIF 6266 ! 6267 !-- Left-right exchange if more than one PE subdomain in the x-direction. 6268 !-- Note that in case of 3-D 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,icrw-1), 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 6261 IF ( bc_dirichlet_l ) THEN 6262 workarrc_sn(0:cg%nz+1,0:2,iclw:icrw-1) & 6263 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw:icrw-1) 6264 ELSE IF ( bc_dirichlet_r ) THEN 6265 workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw) & 6266 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw) 6285 6267 ELSE 6286 6268 workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw-1) & 6287 6269 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1) 6288 6270 ENDIF 6289 6290 IF ( var == 'u' ) THEN 6271 ! 6272 !-- Left-right exchange if more than one PE subdomain in the x-direction. 6273 !-- Note that in case of 3-D nesting the left (pleft == MPI_PROC_NULL) and 6274 !-- right (pright == MPI_PROC_NULL) boundaries are not exchanged because 6275 !-- the nest domain is not cyclic. 6276 !-- From left to right 6277 CALL MPI_SENDRECV( workarrc_sn(0,0,iclw+1), 1, & 6278 workarrc_sn_exchange_type, pleft, 0, & 6279 workarrc_sn(0,0,icrw), 1, & 6280 workarrc_sn_exchange_type, pright, 0, & 6281 comm2d, status, ierr ) 6282 ! 6283 !-- From right to left 6284 CALL MPI_SENDRECV( workarrc_sn(0,0,icrw-1), 1, & 6285 workarrc_sn_exchange_type, pright, 1, & 6286 workarrc_sn(0,0,iclw), 1, & 6287 workarrc_sn_exchange_type, pleft, 1, & 6288 comm2d, status, ierr ) 6289 #endif 6290 ELSE 6291 workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw-1) & 6292 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1) 6293 ENDIF 6294 6295 IF ( var == 'u' ) THEN 6291 6296 6292 DO l = iclw+1, icrw 6293 DO n = 0, kct 6294 ! 6295 !-- First substitute only the matching-node values 6297 DO l = iclw, icrw-1 6298 DO n = 0, kct 6299 ! 6300 !-- Use averages of the neighbouring matching grid-line values 6301 DO i = ifl(l), ifl(l+1) 6302 f(kfl(n):kfu(n),jbc,i) = 0.5_wp * ( workarrc_sn(n,mw,l) & 6303 + workarrc_sn(n,mw,l+1) ) 6304 ENDDO 6305 ! 6306 !-- Then set the values along the matching grid-lines 6307 IF ( MOD( ifl(l), igsr ) == 0 ) THEN 6296 6308 f(kfl(n):kfu(n),jbc,ifl(l)) = workarrc_sn(n,mw,l) 6297 6309 ENDIF 6310 6311 ENDDO 6312 ENDDO 6313 ! 6314 !-- Finally, set the values along the last matching grid-line 6315 IF ( MOD( ifl(icrw), igsr ) == 0 ) THEN 6316 DO n = 0, kct 6317 f(kfl(n):kfu(n),jbc,ifl(icrw)) = workarrc_sn(n,mw,icrw) 6318 ENDDO 6319 ENDIF 6320 ! 6321 !-- A gap may still remain in some cases if the subdomain size is not 6322 !-- divisible by the grid-spacing ratio. In such a case, fill the 6323 !-- gap. Note however, this operation may produce some additional 6324 !-- momentum conservation error. 6325 IF ( ifl(icrw) < nxr ) THEN 6326 DO n = 0, kct 6327 DO i = ifl(icrw)+1, nxr 6328 f(kfl(n):kfu(n),jbc,i) = f(kfl(n):kfu(n),jbc,ifl(icrw)) 6298 6329 ENDDO 6299 6330 ENDDO 6300 6301 DO l = iclw+1, icrw-1 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)) ) 6308 ENDDO 6309 6331 ENDIF 6332 6333 ELSE IF ( var == 'w' ) THEN 6334 6335 DO l = iclw, icrw 6336 DO n = 0, kct + 1 ! It is important to go up to kct+1 6337 ! 6338 !-- First substitute only the matching-node values 6339 f(kfu(n),jbc,ifl(l):ifu(l)) = workarrc_sn(n,mw,l) 6340 6341 ENDDO 6342 ENDDO 6343 6344 DO l = iclw, icrw 6345 DO n = 1, kct + 1 ! It is important to go up to kct+1 6346 ! 6347 !-- Then fill up the nodes in between with the averages 6348 DO k = kfu(n-1)+1, kfu(n)-1 6349 f(k,jbc,ifl(l):ifu(l)) = 0.5_wp * ( f(kfu(n-1),jbc,ifl(l):ifu(l)) & 6350 + f(kfu(n),jbc,ifl(l):ifu(l)) ) 6310 6351 ENDDO 6352 6311 6353 ENDDO 6312 6313 ELSE IF ( var == 'w' ) THEN 6314 6315 DO l = iclw+1, icrw-1 6316 DO n = 0, kct + 1 ! It is important to go up to kct+1 6317 ! 6318 !-- First substitute only the matching-node 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, icrw-1 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(n-1)+1, kfu(n)-1 6329 f(k,jbc,ifl(l):ifu(l)) = 0.5_wp * ( f(kfu(n-1),jbc,ifl(l):ifu(l)) & 6330 + f(kfu(n),jbc,ifl(l):ifu(l)) ) 6331 ENDDO 6332 6333 ENDDO 6334 ENDDO 6335 6336 ENDIF 6354 ENDDO 6337 6355 6338 6356 ELSE ! v or scalars 6339 6357 6340 DO l = iclw +1, icrw-16358 DO l = iclw, icrw 6341 6359 DO n = 0, kct 6342 6360 6343 6361 DO i = ifl(l), ifu(l) 6344 6362 DO k = kfl(n), kfu(n) 6345 f(k,jbc,i) = fc(n,m,l)6363 f(k,jbc,i) = workarrc_sn(n,mw,l) 6346 6364 ENDDO 6347 6365 ENDDO … … 6374 6392 !-- Interpolation of ghost-node values used as the child-domain boundary 6375 6393 !-- conditions. This subroutine handles the top boundary. 6376 !-- This subroutine is based on trilinear interpolation.6377 6378 6394 IMPLICIT NONE 6379 6395 … … 6459 6475 REAL(wp) :: rcorr_ijk !< 6460 6476 6461 real(wp), parameter :: c1 = 2.0_wp / 6.0_wp6462 real(wp), parameter :: c2 = 5.0_wp / 6.0_wp6463 real(wp), parameter :: c3 = -1.0_wp / 6.0_wp6477 ! real(wp), parameter :: c1 = 2.0_wp / 6.0_wp 6478 ! real(wp), parameter :: c2 = 5.0_wp / 6.0_wp 6479 ! real(wp), parameter :: c3 = -1.0_wp / 6.0_wp 6464 6480 6465 6481 … … 6481 6497 var_flag = 0 6482 6498 ENDIF 6483 n = kc(k) + noff 6499 n = kc(k) + noff ! Chance this to get rid of kc. 6484 6500 nw = noff 6485 6501 ! write(9,"(a1,2x,5(i3,2x))") var, k, kc(k), noff, n, nw … … 6551 6567 #endif 6552 6568 6553 IF ( var == 'w' ) THEN 6554 ! DO l = iclw+1, icrw-1 6555 ! DO m = jcsw+1, jcnw-1 6556 DO l = iclw, icrw ! Also fill up beyond edges 6557 DO m = jcsw, jcnw ! Also fill up beyond edges 6558 6569 IF ( var == 'u' ) THEN 6570 6571 DO l = iclw, icrw-1 6572 DO m = jcsw, jcnw 6573 ! 6574 !-- Use averages of the neighbouring matching grid-line values 6575 DO i = ifl(l), ifl(l+1) 6576 f(k,jfl(m):jfu(m),i) = 0.5_wp * ( workarrc_t(nw,m,l) & 6577 + workarrc_t(nw,m,l+1) ) 6578 ENDDO 6579 ! 6580 !-- Then set the values along the matching grid-lines 6581 IF ( MOD( ifl(l), igsr ) == 0 ) THEN 6582 f(k,jfl(m):jfu(m),ifl(l)) = workarrc_t(nw,m,l) 6583 ENDIF 6584 6585 ENDDO 6586 ENDDO 6587 ! 6588 !-- Finally, set the values along the last matching grid-line 6589 IF ( MOD( ifl(icrw), igsr ) == 0 ) THEN 6590 DO m = jcsw, jcnw 6591 f(k,jfl(m):jfu(m),ifl(icrw)) = workarrc_t(nw,m,icrw) 6592 ENDDO 6593 ENDIF 6594 ! 6595 !-- A gap may still remain in some cases if the subdomain size is not 6596 !-- divisible by the grid-spacing ratio. In such a case, fill the 6597 !-- gap. Note however, this operation may produce some additional 6598 !-- momentum conservation error. 6599 IF ( ifl(icrw) < nxr ) THEN 6600 DO m = jcsw, jcnw 6601 DO i = ifl(icrw)+1, nxr 6602 f(k,jfl(m):jfu(m),i) = f(k,jfl(m):jfu(m),ifl(icrw)) 6603 ENDDO 6604 ENDDO 6605 ENDIF 6606 6607 ELSE IF ( var == 'v' ) THEN 6608 6609 DO l = iclw, icrw 6610 DO m = jcsw, jcnw-1 6611 ! 6612 !-- Use averages of the neighbouring matching grid-line values 6613 DO j = jfl(m), jfl(m+1) 6614 f(k,j,ifl(l):ifu(l)) = 0.5_wp * ( workarrc_t(nw,m,l) & 6615 + workarrc_t(nw,m+1,l) ) 6616 ENDDO 6617 ! 6618 !-- Then set the values along the matching grid-lines 6619 IF ( MOD( jfl(m), jgsr ) == 0 ) THEN 6620 f(k,jfl(m),ifl(l):ifu(l)) = workarrc_t(nw,m,l) 6621 ENDIF 6622 6623 ENDDO 6624 6625 ENDDO 6626 ! 6627 !-- Finally, set the values along the last matching grid-line 6628 IF ( MOD( jfl(jcnw), jgsr ) == 0 ) THEN 6629 DO l = iclw, icrw 6630 f(k,jfl(jcnw),ifl(l):ifu(l)) = workarrc_t(nw,jcnw,l) 6631 ENDDO 6632 ENDIF 6633 ! 6634 !-- A gap may still remain in some cases if the subdomain size is not 6635 !-- divisible by the grid-spacing ratio. In such a case, fill the 6636 !-- gap. Note however, this operation may produce some additional 6637 !-- momentum conservation error. 6638 IF ( jfl(jcnw) < nyn ) THEN 6639 DO l = iclw, icrw 6640 DO j = jfl(jcnw)+1, nyn 6641 f(k,j,ifl(l):ifu(l)) = f(k,jfl(jcnw),ifl(l):ifu(l)) 6642 ENDDO 6643 ENDDO 6644 ENDIF 6645 6646 ELSE ! w or any scalar 6647 6648 DO l = iclw, icrw 6649 DO m = jcsw, jcnw 6650 ! 6651 !-- 3rd-order upwind biased interpolation. 6652 ! IF ( w(k,jfl(m),ifl(l)) < 0.0_wp ) THEN 6653 ! f_interp = c1 * workarrc_t(0,m,l) + c2 * workarrc_t(1,m,l) + c3 * workarrc_t(2,m,l) 6654 ! ELSE 6655 ! f_interp = workarrc_t(nw,m,l) 6656 ! ENDIF 6559 6657 DO i = ifl(l), ifu(l) 6560 6658 DO j = jfl(m), jfu(m) … … 6565 6663 ENDDO 6566 6664 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, jcnw-1 6573 ! 6574 !-- 1st-order upwind interpolation 6575 f_interp = workarrc_t(nw,m,l) 6576 ! 6577 !-- 3rd-order 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 matching-node 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 6599 DO l = iclw+1, icrw-1 6600 DO m = jcsw+1, jcnw-1 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, icrw-1 6615 DO m = jcsw+1, jcnw ! Note that we have to go up to jcnw here. 6616 ! 6617 !-- 1st-order upwind interpolation 6618 f_interp = workarrc_t(nw,m,l) 6619 ! 6620 !-- 3rd-order 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 matching-node 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, icrw-1 6642 DO m = jcsw+1, jcnw-1 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, icrw-1 6656 DO m = jcsw+1, jcnw-1 6657 ! 6658 !-- 1st-order upwind interpolation 6659 f_interp = workarrc_t(nw,m,l) 6660 ! 6661 !-- 3rd-order 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 6669 DO i = ifl(l), ifu(l) 6670 DO j = jfl(m), jfu(m) 6671 f(k,j,i) = f_interp 6672 ENDDO 6673 ENDDO 6674 6675 ENDDO 6676 ENDDO 6677 ENDIF 6665 6666 ENDIF ! var 6678 6667 ! 6679 6668 !-- Just fill up the redundant second ghost-node layer in case of var == w.
Note: See TracChangeset
for help on using the changeset viewer.
