source: palm/trunk/SOURCE/init_rankine.f90 @ 9

 SUBROUTINE init_rankine

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: init_rankine.f90 4 2007-02-13 11:33:16Z raasch $
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.11  2005/03/26 20:38:49  raasch
! Arguments for non-cyclic boundary conditions added to argument list of
! routine exchange_horiz
!
! Revision 1.1  1997/08/11 06:18:43  raasch
! Initial revision
!
!
! Description:
! ------------
! Initialize a (nondivergent) Rankine eddy with a vertical axis in order to test 
! the advection terms and the pressure solver.
!------------------------------------------------------------------------------!

    USE arrays_3d
    USE constants
    USE grid_variables
    USE indices
    USE control_parameters

    IMPLICIT NONE

    INTEGER ::  i, ic, j, jc, k, kc1, kc2
    REAL    ::  alpha, betrag, radius, rc, uw, vw, x, y

!
!-- Default: eddy radius rc, eddy strength z,
!--          position of eddy centre: ic, jc, kc1, kc2
    rc  =  4.0 * dx
    ic  =  ( nx+1 ) / 2
    jc  =  ic
    kc1 = nzb
    kc2 = nzt+1

!
!-- Compute the u-component.
    DO  i = nxl, nxr
       DO  j = nys, nyn
          x = ( i - ic - 0.5 ) * dx
          y = ( j - jc ) * dy
          radius = SQRT( x**2 + y**2 )
          IF ( radius <= 2.0 * rc )  THEN
             betrag = radius / ( 2.0 * rc ) * 0.08
          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
          ELSE
             betrag = 0.0
          ENDIF

          IF ( x == 0.0 )  THEN
             IF ( y > 0.0 )  THEN
                alpha = pi / 2.0
             ELSEIF ( y < 0.0 )  THEN
                alpha = 3.0 * pi / 2.0
             ENDIF
          ELSE
             IF ( x < 0.0 )  THEN
                alpha = ATAN( y / x ) + pi
             ELSE
                IF ( y < 0.0 )  THEN
                   alpha = ATAN( y / x ) + 2.0 * pi
                ELSE
                   alpha = ATAN( y / x )
                ENDIF
             ENDIF
          ENDIF

          uw = -SIN( alpha ) * betrag

          DO  k = kc1, kc2
             u(k,j,i) = u(k,j,i) + uw
          ENDDO
       ENDDO
    ENDDO

!
!-- Compute the v-component.
    DO  i = nxl, nxr
       DO  j = nys, nyn
          x = ( i - ic ) * dx
          y = ( j - jc - 0.5) * dy
          radius = SQRT( x**2 + y**2 )
          IF ( radius <= 2.0 * rc )  THEN
             betrag = radius / ( 2.0 * rc ) * 0.08
          ELSEIF ( radius > 2.0 * rc  .AND.  radius < 8.0 * rc )  THEN
             betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 )
          ELSE
             betrag = 0.0
          ENDIF

          IF ( x == 0.0 )  THEN
             IF ( y > 0.0 )  THEN
                alpha = pi / 2.0
             ELSEIF ( y < 0.0 )  THEN
                alpha = 3.0 * pi / 2.0
             ENDIF
          ELSE
             IF ( x < 0.0 )  THEN
                alpha = ATAN( y / x ) + pi
             ELSE
                IF ( y < 0.0 )  THEN
                   alpha = ATAN( y / x ) + 2.0 * pi
                ELSE
                   alpha = ATAN( y / x )
                ENDIF
             ENDIF
          ENDIF

          vw = COS( alpha ) * betrag

          DO  k = kc1, kc2
             v(k,j,i) = v(k,j,i) + vw
          ENDDO
       ENDDO
    ENDDO

!
!-- Exchange of boundary values for the velocities.
    CALL exchange_horiz( u, uxrp,    0 )
    CALL exchange_horiz( v,    0, vynp )
!
!-- Make velocity field nondivergent.
    n_sor = nsor_ini
    CALL pres
    n_sor = nsor

 END SUBROUTINE init_rankine