source: palm/trunk/SOURCE/advec_u_ups.f90 @ 803

 SUBROUTINE advec_u_ups 

!------------------------------------------------------------------------------!
! Current revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: advec_u_ups.f90 484 2010-02-05 07:36:54Z franke $
!
! 164 2008-05-15 08:46:15Z raasch
! Arguments removed from transpose routines
!
! February 2007
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.6  2004/04/30 08:03:19  raasch
! Enlarged transposition arrays introduced
!
! Revision 1.1  1999/02/05 08:49:08  raasch
! Initial revision
!
!
! Description:
! ------------
! Upstream-Spline advection of the u velocity-component. The advection process
! is divided into three subsequent steps, one for each of the dimensions. The
! results is stored as a tendency in array tend. The computation of the cubic
! splines and the possible execution of the Long-filter require that all grid
! points of the relevant dimension are available. For model runs on more than
! one PE therefore both the advected and the advecting quantities are
! transposed accordingly.
!
! Internally used arrays:
! v_ad = scalar quantity to be advected, initialised = u at the beginning,
!        also being used as temporary storage after each time step
! d    = advecting component (u, v, or w)
!------------------------------------------------------------------------------!

    USE advection
    USE arrays_3d
    USE cpulog
    USE grid_variables
    USE indices
    USE interfaces
    USE control_parameters

    IMPLICIT NONE

    INTEGER :: i, j, k
    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  v_ad

    
    CALL cpu_log( log_point_s(17), 'advec_u_ups', 'start' )

#if defined( __parallel )

!
!-- Advection of u in x-direction:
!-- Store u in temporary array v_ad (component to be advected, boundaries
!-- are not used because they disturb the transposition)
    ALLOCATE( v_ad(nzb+1:nzta,nys:nyna,nxl:nxra) )
    v_ad = 0.0
    v_ad(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)

!
!-- Enlarge the size of tend, used as a working array for the transpositions
    IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
       DEALLOCATE( tend )
       ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
    ENDIF

!
!-- Transpose the component to be advected: z --> x
    CALL transpose_zx( v_ad, tend, v_ad )

!
!-- Advecting component (d) = component to be advected (v_ad) (variable d is
!-- used for storage, because it is the only one having suitable dimensions).
!-- NOTE: here x is the first dimension and lies completely on the PE.
    d = v_ad - u_gtrans

#else

!
!-- Advection of u in x-direction:
!-- Store u in temporary array v_ad (component to be advected)
    ALLOCATE( v_ad(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
    v_ad(:,:,:) = u(:,:,:)

!
!-- Advecting component (d) = component to be advected (u) (variable d is used
!-- for storage, because it is the only one having suitable dimensions. This is 
!-- done for for reasons of compatibility with the parallel part.)
    d(:,:,:) = u(nzb+1:nzt,nys:nyn,nxl:nxr) - u_gtrans

#endif

!
!-- Upstream-Spline advection of u in x-direction. Array tend comes out
!-- as v_ad before the advection step including cyclic boundaries.
!-- It is needed for the long filter.
    CALL spline_x( v_ad, d, 'u' )

!
!-- Advection of u in y-direction:
!-- advecting component (v) must be averaged out on the u grid
    DO  i = nxl, nxr
       DO  j = nys, nyn
          DO  k = nzb+1, nzt
              d(k,j,i) = 0.25 * ( v(k,j,i-1) + v(k,j+1,i-1) + &
                                  v(k,j,i)   + v(k,j+1,i) ) - v_gtrans
          ENDDO
       ENDDO
    ENDDO

#if defined( __parallel )

!
!-- Transpose the advecting component: z --> y
    CALL transpose_zx( d, tend, d )
    CALL transpose_xy( d, tend, d )

!
!-- Transpose the component to be advected: x --> y
    CALL transpose_xy( v_ad, tend, v_ad )

#endif

!
!-- Upstream-Spline advection of u in y-direction
    CALL spline_y( v_ad, d, 'u' )

!
!-- Advection of u in z-direction:
!-- the advecting component (w) must be averaged out on the u grid
!-- (weighted for non-equidistant grid)
    DO  i = nxl, nxr
       DO  j = nys, nyn
          DO  k = nzb+1, nzt
             d(k,j,i) =  ( 0.5 * ( w(k,j,i)   + w(k,j,i-1)   ) * &
                                 ( zu(k) - zw(k-1) )             &
                         + 0.5 * ( w(k-1,j,i) + w(k-1,j,i-1) ) * &
                                 ( zw(k) - zu(k)   )             &
                         ) * ddzw(k)
          ENDDO
       ENDDO
    ENDDO

#if defined( __parallel )

!
!-- Transpose the component to be advected: y --> z (= y --> x + x --> z)
    CALL transpose_yx( v_ad, tend, v_ad )
    CALL transpose_xz( v_ad, tend, v_ad )

!
!-- Resize tend to its normal size
    IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
       DEALLOCATE( tend )
       ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
    ENDIF

#endif

!
!-- Upstream-Spline advection of u in z-direction
    CALL spline_z( v_ad, d, dzu, spl_tri_zu, 'u' )

!
!-- Compute the tendency term
    DO  i = nxl, nxr
       DO  j = nys, nyn
          DO  k = nzb+1, nzt
             tend(k,j,i) = ( v_ad(k,j,i) - u(k,j,i) ) / dt_3d 
          ENDDO
       ENDDO
    ENDDO

    DEALLOCATE( v_ad )

    CALL cpu_log( log_point_s(17), 'advec_u_ups', 'stop' )

 END SUBROUTINE advec_u_ups