Ignore:
Timestamp:
Jan 30, 2019 12:58:13 PM (3 years ago)
Author:
hellstea
Message:

Checks for parent / child grid line matching introduced

File:
1 edited

Legend:

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

    r3697 r3708  
    2525! -----------------
    2626! $Id$
     27! Checks for parent / child grid line matching introduced.
     28! Interpolation of nest-boundary-tangential velocity components revised.
     29!
     30! 3697 2019-01-24 17:16:13Z hellstea
    2731! Bugfix: upper k-bound in the child initialization interpolation
    2832! pmci_interp_1sto_all corrected.
     
    14751479!--    included in the interpolation algorithms.
    14761480       CALL pmci_init_anterp_tophat
     1481!
     1482!--    Check that the child and parent grid lines match
     1483       CALL pmci_check_grid_matching
    14771484
    14781485    ENDIF
     
    19281935            CEILING( parentdzmax / dzmin )
    19291936       ALLOCATE( celltmpd(1:acsize) )
    1930      
     1937
    19311938    END SUBROUTINE pmci_init_interp_tril
    19321939   
     
    19411948       igsr = NINT( cg%dx / dx, iwp )
    19421949       jgsr = NINT( cg%dy / dy, iwp )
    1943        kgsr = NINT( cg%dzu(cg%nz+1) / dzu(nzt+1), iwp )
     1950       kgsr = NINT( cg%dzw(1) / dzw(1), iwp )
    19441951       write(9,"('igsr, jgsr, kgsr: ',3(i3,2x))") igsr, jgsr, kgsr
    19451952       flush(9)
     
    38993906    END SUBROUTINE pmci_init_anterp_tophat
    39003907
    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.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
     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 grid-spacing 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 grid-spacing 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 grid-spacing 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 grid-spacing ratio ',           &
     3973               '( parent dz / child dz ) must have an integer value for each z-level'
     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 
    39023981 END SUBROUTINE pmci_setup_child
    39033982
     
    45574636             DO  m = mb, me
    45584637                DO n = 0, kct + 1 
    4559                
     4638
    45604639                   DO i = ifl(l), ifu(l)
    45614640                      DO  j = jfl(m), jfl(m+1)-1
     
    45944673             DO  m = mb, me
    45954674                DO n = 0, kct + 1
    4596                    
     4675                                      
    45974676                   DO i = ifl(l), ifu(l)
    4598                       DO  j = jfl(m), jfu(m)
     4677                      DO  j = jfl(m), jfu(m)                         
    45994678                         DO  k = kfl(n), MIN( kfu(n), nzt+1 )
    46004679                            f(k,j,i) = fc(n,m,l)
     
    57855864      REAL(wp) ::  rcorr        !< Average reversibility correction for the whole anterpolation cell
    57865865      REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid 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
    57875870
    57885871!
     
    58585941      ENDIF
    58595942
    5860       IF  ( var == 'u' )  THEN
    5861 !AH!
    5862 !AH!--      Substitute the necessary parent-grid 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:jcnw-1,0:2)                           &
    5868 !AH                    = fc(0:cg%nz+1,jcsw:jcnw-1,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:jcnw-1,0:2)                         &
    5874 !AH                    = fc(0:cg%nz+1,jcsw+1:jcnw-1,lbeg:lbeg+2)
    5875 !AH            ENDIF
    5876 !AH!
    5877 !AH!--         South-north exchange if more than one PE subdomain in the y-direction.
    5878 !AH!--         Note that in case of 3-D 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,jcnw-1,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, jcnw-1
    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 parent-grid 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: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
     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 matching-node 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, 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) )
    59066001                  ENDDO
     6002                 
    59076003               ENDDO
    5908 
    59096004            ENDDO
    5910          ENDDO
    5911 
    5912       ELSE IF  ( var == 'v' )  THEN
    5913          
    5914          DO  m = jcsw+1, jcnw-1
    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, 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) )
    59206024                  ENDDO
     6025                 
    59216026               ENDDO
    5922 
    59236027            ENDDO
    5924          ENDDO
    5925 
    5926       ELSE IF  ( var == 'w' )  THEN
    5927 
    5928          DO  m = jcsw+1, jcnw-1
    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 n-loop starts from n=1 instead of 0
    5933                   DO  k = kfu(n-1)+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
    59426032         
    59436033         DO  m = jcsw+1, jcnw-1
     
    60816171      REAL(wp) ::  rcorr        !< Average reversibility correction for the whole anterpolation cell
    60826172      REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid 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
    60836177
    60846178!
     
    61546248      ENDIF
    61556249
    6156       IF  ( var == 'v' )  THEN
    6157 !AH!
    6158 !AH!--      Substitute the necessary parent-grid 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:icrw-1)                           &
    6164 !AH                    = fc(0:cg%nz+1,mbeg:mbeg+2,iclw:icrw-1)
    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:icrw-1)                         &
    6170 !AH                    = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1)
    6171 !AH            ENDIF
    6172 !AH!
    6173 !AH!--         Left-right exchange if more than one PE subdomain in the x-direction.
    6174 !AH!--         Note that in case of 3-D 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,icrw-1), 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:icrw-1)                            &
    6193 !AH                 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1)
    6194 !AH         ENDIF
    6195 
    6196          DO  l = iclw+1, icrw-1
    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 parent-grid 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: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     
     6285         ELSE
     6286            workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw-1)                            &
     6287                 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1)
     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 matching-node 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, 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)) )
    62026308                  ENDDO
     6309
    62036310               ENDDO
    6204 
    62056311            ENDDO
    6206          ENDDO
    6207 
    6208       ELSE IF  ( var == 'u' )  THEN
    6209          
    6210          DO  l = iclw+1, icrw-1
    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, 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)) )
    62206331                  ENDDO
     6332
    62216333               ENDDO
    6222 
    62236334            ENDDO
    6224          ENDDO
    6225 
    6226       ELSE IF  ( var == 'w' )  THEN
    6227 
    6228          DO  l = iclw+1, icrw-1
    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 n-loop starts from n=1 instead of 0
    6233                   DO  k = kfu(n-1)+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
    62426339         
    62436340         DO  l = iclw+1, icrw-1
     
    63586455      REAL(wp) ::  fkpjp       !<
    63596456      REAL(wp) ::  fkp         !<
     6457      REAL(wp) ::  f_interp    !<
    63606458      REAL(wp) ::  rcorr       !<
    63616459      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
    63626464
    63636465
     
    63806482      ENDIF
    63816483      n  = kc(k) + noff
    6382       nw = 1
    6383 
    6384       IF  ( var == 'w' )  THEN
    6385 !AH!
    6386 !AH!--      Substitute the necessary parent-grid 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?w-index 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!--      Left-right exchange if more than one PE subdomain in the x-direction.
    6411 !AH!--      Note that in case of 3-D 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,icrw-1), 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!--      South-north exchange if more than one PE subdomain in the y-direction.
    6432 !AH!--      Note that in case of 3-D 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,jcnw-1,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 parent-grid 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?w-index 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!--   Left-right exchange if more than one PE subdomain in the x-direction.
     6512!--   Note that in case of 3-D 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,icrw-1), 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!--   South-north exchange if more than one PE subdomain in the y-direction.
     6533!--   Note that in case of 3-D 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,jcnw-1,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, 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
     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, 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
    64536599         DO  l = iclw+1, icrw-1
    64546600            DO  m = jcsw+1, jcnw-1
    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, 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             
    64566669               DO  i = ifl(l), ifu(l)
    64576670                  DO  j = jfl(m), jfu(m)
    6458                      f(k,j,i) = fc(n,m,l)
     6671                     f(k,j,i) = f_interp
    64596672                  ENDDO
    64606673               ENDDO
     
    64626675            ENDDO
    64636676         ENDDO
    6464 
    6465       ELSE IF  ( var == 'u' )  THEN
    6466          
    6467          DO  l = iclw+1, icrw-1
    6468             DO  m = jcsw+1, jcnw-1
    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, icrw-1
    6482             DO  m = jcsw+1, jcnw-1
    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, icrw-1
    6496             DO  m = jcsw+1, jcnw-1
    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
    65086678!
    65096679!--   Just fill up the redundant second ghost-node layer in case of var == w.
Note: See TracChangeset for help on using the changeset viewer.