Changeset 3741


Ignore:
Timestamp:
Feb 13, 2019 4:24:49 PM (2 years ago)
Author:
hellstea
Message:

Bugfix in nest initialization and interpolations

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/pmc_interface_mod.f90

    r3708 r3741  
    2525! -----------------
    2626! $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
    2733! Checks for parent / child grid line matching introduced.
    2834! Interpolation of nest-boundary-tangential velocity components revised.
     
    15041510       INTEGER(iwp) :: jauxs    !<
    15051511       INTEGER(iwp) :: jauxn    !<
    1506        REAL(wp) ::  loffset     !<
    1507        REAL(wp) ::  noffset     !<
    1508        REAL(wp) ::  roffset     !<
    1509        REAL(wp) ::  soffset     !<
    15101512       REAL(wp) ::  xexl        !< Parent-grid array exceedance behind the left edge of the child PE subdomain
    15111513       REAL(wp) ::  xexr        !< Parent-grid array exceedance behind the right edge of the child PE subdomain
     
    15801582!--    Left
    15811583       IF  ( bc_dirichlet_l )  THEN
    1582           loffset = MOD( coord_x(nxl), cg%dx )
    1583           xexl  = 2 * cg%dx + loffset
     1584          xexl  = 2 * cg%dx
    15841585          iauxl = 0
    15851586       ELSE
     
    15971598!--    Right
    15981599       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
    16011601          iauxr = 0 
    16021602       ELSE
     
    16071607       DO  i = cg%nx, 0 , -1
    16081608          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 ) )
    16101610             EXIT
    16111611          ENDIF
     
    16141614!--    South
    16151615       IF  ( bc_dirichlet_s )  THEN
    1616           soffset = MOD( coord_y(nys), cg%dy )
    1617           yexs  = 2 * cg%dy + soffset
     1616          yexs  = 2 * cg%dy
    16181617          jauxs = 0 
    16191618       ELSE
     
    16311630!--    North
    16321631       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
    16351633          jauxn = 0
    16361634       ELSE
     
    16411639       DO  j = cg%ny, 0 , -1
    16421640          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 ) )
    16441642             EXIT
    16451643          ENDIF
     
    35563554!
    35573555!--    i-indices of u for each ii-index value
    3558 !--    ii=icr is redundant for anterpolation
    3559        tolerance = 0.000001_wp * dx
    35603556       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.
    35693561          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 )
    35713563             i = i + 1
    35723564          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)
    35973568!
    35983569!--       Print out the index bounds for checking and debugging purposes
     
    36053576!
    36063577!--    i-indices of others for each ii-index value
    3607 !--    ii=icr is redundant for anterpolation
    36083578       istart = nxlg
    3609 !AH       DO  ii = icl, icr
    3610 !AH       DO  ii = icl-1, icr+1
    36113579       DO  ii = icla, icra
    36123580          i = istart
     
    36333601!
    36343602!--    j-indices of v for each jj-index value
    3635 !--    jj=jcn is redundant for anterpolation
    3636        tolerance = 0.000001_wp * dy
    36373603       jstart = nysg
    3638 !AH       DO  jj = jcs, jcn
    3639 !AH       DO  jj = jcs-1, jcn+1
    36403604       DO  jj = jcsa, jcna
    36413605!
    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.
    36453608          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 )
    36473610             j = j + 1
    36483611          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)
    36733615!
    36743616!--       Print out the index bounds for checking and debugging purposes
     
    36803622!
    36813623!--    j-indices of others for each jj-index value
    3682 !--    jj=jcn is redundant for anterpolation
    36833624       jstart = nysg
    3684 !AH       DO  jj = jcs, jcn
    3685 !AH       DO  jj = jcs-1, jcn+1
    36863625       DO  jj = jcsa, jcna
    36873626          j = jstart
     
    37093648!--    k-indices of w for each kk-index value
    37103649!--    Note that anterpolation index limits are needed also for the top boundary
    3711 !--    ghost cell level because of the reversibility correction in the interpolation.
     3650!--    ghost cell level because they are used also in the interpolation.
    37123651       kstart  = 0
    37133652       kflw(0) = 0
    37143653       kfuw(0) = 0
    3715        tolerance = 0.000001_wp * dzw(1)
    37163654       DO kk = 1, cg%nz+1
    37173655!
    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.
    37213658          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 ) )
    37233660             k = k + 1
    37243661          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
    37513667          write(9,"('pmci_init_anterp_tophat, kk, kflw, kfuw: ', 4(i4,2x), 2(e12.5,2x))") &
    37523668               kk, kflw(kk), kfuw(kk), nzt,  cg%zu(kk), cg%zw(kk)
     
    37613677!
    37623678!--    Note that anterpolation index limits are needed also for the top boundary
    3763 !--    ghost cell level because of the reversibility correction in the interpolation.
     3679!--    ghost cell level because they are used also in the interpolation.
    37643680       DO  kk = 1, cg%nz+1
    37653681          k = kstart
     
    37823698       kflo(cg%nz+1) = nzt+kgsr ! Obsolete
    37833699       kfuo(cg%nz+1) = nzt+kgsr ! Obsolete
    3784 
     3700!
     3701!--       Print out the index bounds for checking and debugging purposes
    37853702       DO  kk = 1, cg%nz+1
    37863703          write(9,"('pmci_init_anterp_tophat, kk, kflo, kfuo: ', 4(i4,2x), 2(e12.5,2x))") &
     
    37893706       ENDDO
    37903707       write(9,*)
    3791 !AH
    3792 
    37933708!
    37943709!--    Precomputation of number of fine-grid nodes inside parent-grid cells.
     
    39103825    SUBROUTINE pmci_check_grid_matching
    39113826!
    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.
    39133830       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
    39223843
    39233844
     
    39263847       non_int_gsr_y = 0
    39273848       non_int_gsr_z = 0
     3849       too_narrow_pesd_x = 0
     3850       too_narrow_pesd_y = 0
    39283851
    39293852       IF  ( myid == 0 )  THEN
     
    39443867          ENDDO
    39453868
     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                         
    39463874       ENDIF
    39473875
     
    39733901               '( parent dz / child dz ) must have an integer value for each z-level'
    39743902          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       
    39763920
    39773921     END SUBROUTINE pmci_check_grid_matching
     
    45214465!
    45224466!--    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/t
    45254467       IMPLICIT NONE
    45264468
     
    45584500
    45594501       INTEGER(iwp) ::  i        !<
     4502       INTEGER(iwp) ::  ib       !<
     4503       INTEGER(iwp) ::  ie       !<
     4504       INTEGER(iwp) ::  ifirst   !<
     4505       INTEGER(iwp) ::  ilast    !<
    45604506       INTEGER(iwp) ::  j        !<
     4507       INTEGER(iwp) ::  jb       !<
     4508       INTEGER(iwp) ::  je       !<
     4509       INTEGER(iwp) ::  jfirst   !<
     4510       INTEGER(iwp) ::  jlast    !<
    45614511       INTEGER(iwp) ::  k        !<
    45624512       INTEGER(iwp) ::  l        !<
     
    45744524       mb = jcs
    45754525       me = jcn
     4526       ifirst = nxl
     4527       ilast  = nxr
     4528       jfirst = nys
     4529       jlast  = nyn
    45764530
    45774531       IF ( nesting_mode /= 'vertical' )  THEN
    45784532          IF ( bc_dirichlet_l )  THEN
    4579              lb = icl + 1
     4533             lb     = icl + 1
     4534             ifirst = nxl - 1
    45804535!
    45814536!--          For u, nxl is a ghost node, but not for the other variables
    45824537             IF ( var == 'u' )  THEN
    4583                 lb = icl + 2
     4538                lb     = icl + 2
     4539                ifirst = nxl
    45844540             ENDIF
    45854541          ENDIF
    45864542          IF ( bc_dirichlet_s )  THEN
    4587              mb = jcs + 1
     4543             mb     = jcs + 1
     4544             jfirst = nys - 1
    45884545!
    45894546!--          For v, nys is a ghost node, but not for the other variables
    45904547             IF ( var == 'v' )  THEN
    4591                 mb = jcs + 2
     4548                mb     = jcs + 2
     4549                jfirst = nys
    45924550             ENDIF
    45934551          ENDIF
    45944552          IF ( bc_dirichlet_r )  THEN
    4595              le = icr - 1
     4553             le    = icr - 1
     4554             ilast = nxr + 1
    45964555          ENDIF
    45974556          IF ( bc_dirichlet_n )  THEN
    4598              me = jcn - 1
     4557             me    = jcn - 1
     4558             jlast = nyn + 1
    45994559          ENDIF
    46004560       ENDIF
     
    46154575       IF  ( var == 'u' )  THEN
    46164576
     4577          ib = ifl(lb)
     4578          ie = ifl(le+1) - 1
     4579          jb = jfl(mb)
     4580          je = jfu(me)
    46174581          DO  l = lb, le
    46184582             DO  m = mb, me
     
    46334597       ELSE IF  ( var == 'v' )  THEN
    46344598
     4599          ib = ifl(lb)
     4600          ie = ifu(le)
     4601          jb = jfl(mb)
     4602          je = jfl(me+1) - 1
    46354603          DO  l = lb, le
    46364604             DO  m = mb, me
     
    46514619       ELSE IF  ( var == 'w' )  THEN
    46524620
     4621          ib = ifl(lb)
     4622          ie = ifu(le)
     4623          jb = jfl(mb)
     4624          je = jfu(me)
    46534625          DO  l = lb, le
    46544626             DO  m = mb, me
     
    46704642       ELSE   ! scalars
    46714643
     4644          ib = ifl(lb)
     4645          ie = ifu(le)
     4646          jb = jfl(mb)
     4647          je = jfu(me)
    46724648          DO  l = lb, le
    46734649             DO  m = mb, me
     
    46874663
    46884664       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
    46894690
    46904691    END SUBROUTINE pmci_interp_1sto_all
     
    57545755!
    57555756!--   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.
    57595758      IMPLICIT NONE
    57605759
     
    58655864      REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid nodes
    58665865
    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
     5866!      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
    58705869
    58715870!
     
    59405939         var_flag = 0
    59415940      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
    59485945#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)
    59785952         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)
    59815955         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
    59845981         
    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)
    59916013               ENDDO
    59926014            ENDDO
    5993 
    5994             DO  m = jcsw+1, jcnw-1
    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) )
    6001                   ENDDO
     6015         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)
    60026024                 
     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) )
    60036035               ENDDO
     6036                 
    60046037            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
    60306039
    60316040      ELSE   ! u or scalars
    60326041         
    6033          DO  m = jcsw+1, jcnw-1
     6042         DO  m = jcsw, jcnw
    60346043            DO n = 0, kct
    60356044               
    60366045               DO  j = jfl(m), jfu(m)
    60376046                  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)
    60396048                  ENDDO
    60406049               ENDDO
     
    60696078!--   Interpolation of ghost-node values used as the child-domain boundary
    60706079!--   conditions. This subroutine handles the south and north boundaries.
    6071 !--   This subroutine is based on trilinear interpolation.
    6072 
    60736080      IMPLICIT NONE
    60746081
     
    61726179      REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid nodes
    61736180
    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
     6181!      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
    61776184
    61786185!
     
    62476254         var_flag = 0
    62486255      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
    62556260#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)
    62856267         ELSE
    62866268            workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw-1)                            &
    62876269                 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1)
    62886270         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
    62916296         
    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
    62966308                  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))
    62986329               ENDDO
    62996330            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)) )
    63106351               ENDDO
     6352
    63116353            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
    63376355
    63386356      ELSE   ! v or scalars
    63396357         
    6340          DO  l = iclw+1, icrw-1
     6358         DO  l = iclw, icrw
    63416359            DO n = 0, kct
    63426360               
    63436361               DO  i = ifl(l), ifu(l)
    63446362                  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)
    63466364                  ENDDO
    63476365               ENDDO
     
    63746392!--   Interpolation of ghost-node values used as the child-domain boundary
    63756393!--   conditions. This subroutine handles the top boundary.
    6376 !--   This subroutine is based on trilinear interpolation.
    6377 
    63786394      IMPLICIT NONE
    63796395
     
    64596475      REAL(wp) ::  rcorr_ijk   !<
    64606476
    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
     6477!      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
    64646480
    64656481
     
    64816497         var_flag = 0
    64826498      ENDIF
    6483       n  = kc(k) + noff
     6499      n  = kc(k) + noff     ! Chance this to get rid of kc.
    64846500      nw = noff
    64856501!      write(9,"(a1,2x,5(i3,2x))") var, k, kc(k), noff, n, nw
     
    65516567#endif     
    65526568
    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             
    65596657               DO  i = ifl(l), ifu(l)
    65606658                  DO  j = jfl(m), jfu(m)
     
    65656663            ENDDO
    65666664         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
    66786667!
    66796668!--   Just fill up the redundant second ghost-node layer in case of var == w.
Note: See TracChangeset for help on using the changeset viewer.