| [1] | 1 | SUBROUTINE init_rankine |
|---|
| 2 | |
|---|
| 3 | !------------------------------------------------------------------------------! |
|---|
| 4 | ! Actual revisions: |
|---|
| 5 | ! ----------------- |
|---|
| 6 | ! |
|---|
| 7 | ! |
|---|
| 8 | ! Former revisions: |
|---|
| 9 | ! ----------------- |
|---|
| [3] | 10 | ! $Id: init_rankine.f90 4 2007-02-13 11:33:16Z raasch $ |
|---|
| 11 | ! RCS Log replace by Id keyword, revision history cleaned up |
|---|
| 12 | ! |
|---|
| [1] | 13 | ! Revision 1.11 2005/03/26 20:38:49 raasch |
|---|
| 14 | ! Arguments for non-cyclic boundary conditions added to argument list of |
|---|
| 15 | ! routine exchange_horiz |
|---|
| 16 | ! |
|---|
| 17 | ! Revision 1.1 1997/08/11 06:18:43 raasch |
|---|
| 18 | ! Initial revision |
|---|
| 19 | ! |
|---|
| 20 | ! |
|---|
| 21 | ! Description: |
|---|
| 22 | ! ------------ |
|---|
| 23 | ! Initialize a (nondivergent) Rankine eddy with a vertical axis in order to test |
|---|
| 24 | ! the advection terms and the pressure solver. |
|---|
| 25 | !------------------------------------------------------------------------------! |
|---|
| 26 | |
|---|
| 27 | USE arrays_3d |
|---|
| 28 | USE constants |
|---|
| 29 | USE grid_variables |
|---|
| 30 | USE indices |
|---|
| 31 | USE control_parameters |
|---|
| 32 | |
|---|
| 33 | IMPLICIT NONE |
|---|
| 34 | |
|---|
| 35 | INTEGER :: i, ic, j, jc, k, kc1, kc2 |
|---|
| 36 | REAL :: alpha, betrag, radius, rc, uw, vw, x, y |
|---|
| 37 | |
|---|
| 38 | ! |
|---|
| 39 | !-- Default: eddy radius rc, eddy strength z, |
|---|
| 40 | !-- position of eddy centre: ic, jc, kc1, kc2 |
|---|
| 41 | rc = 4.0 * dx |
|---|
| 42 | ic = ( nx+1 ) / 2 |
|---|
| 43 | jc = ic |
|---|
| 44 | kc1 = nzb |
|---|
| 45 | kc2 = nzt+1 |
|---|
| 46 | |
|---|
| 47 | ! |
|---|
| 48 | !-- Compute the u-component. |
|---|
| 49 | DO i = nxl, nxr |
|---|
| 50 | DO j = nys, nyn |
|---|
| 51 | x = ( i - ic - 0.5 ) * dx |
|---|
| 52 | y = ( j - jc ) * dy |
|---|
| 53 | radius = SQRT( x**2 + y**2 ) |
|---|
| 54 | IF ( radius <= 2.0 * rc ) THEN |
|---|
| 55 | betrag = radius / ( 2.0 * rc ) * 0.08 |
|---|
| 56 | ELSEIF ( radius > 2.0 * rc .AND. radius < 8.0 * rc ) THEN |
|---|
| 57 | betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 ) |
|---|
| 58 | ELSE |
|---|
| 59 | betrag = 0.0 |
|---|
| 60 | ENDIF |
|---|
| 61 | |
|---|
| 62 | IF ( x == 0.0 ) THEN |
|---|
| 63 | IF ( y > 0.0 ) THEN |
|---|
| 64 | alpha = pi / 2.0 |
|---|
| 65 | ELSEIF ( y < 0.0 ) THEN |
|---|
| 66 | alpha = 3.0 * pi / 2.0 |
|---|
| 67 | ENDIF |
|---|
| 68 | ELSE |
|---|
| 69 | IF ( x < 0.0 ) THEN |
|---|
| 70 | alpha = ATAN( y / x ) + pi |
|---|
| 71 | ELSE |
|---|
| 72 | IF ( y < 0.0 ) THEN |
|---|
| 73 | alpha = ATAN( y / x ) + 2.0 * pi |
|---|
| 74 | ELSE |
|---|
| 75 | alpha = ATAN( y / x ) |
|---|
| 76 | ENDIF |
|---|
| 77 | ENDIF |
|---|
| 78 | ENDIF |
|---|
| 79 | |
|---|
| 80 | uw = -SIN( alpha ) * betrag |
|---|
| 81 | |
|---|
| 82 | DO k = kc1, kc2 |
|---|
| 83 | u(k,j,i) = u(k,j,i) + uw |
|---|
| 84 | ENDDO |
|---|
| 85 | ENDDO |
|---|
| 86 | ENDDO |
|---|
| 87 | |
|---|
| 88 | ! |
|---|
| 89 | !-- Compute the v-component. |
|---|
| 90 | DO i = nxl, nxr |
|---|
| 91 | DO j = nys, nyn |
|---|
| 92 | x = ( i - ic ) * dx |
|---|
| 93 | y = ( j - jc - 0.5) * dy |
|---|
| 94 | radius = SQRT( x**2 + y**2 ) |
|---|
| 95 | IF ( radius <= 2.0 * rc ) THEN |
|---|
| 96 | betrag = radius / ( 2.0 * rc ) * 0.08 |
|---|
| 97 | ELSEIF ( radius > 2.0 * rc .AND. radius < 8.0 * rc ) THEN |
|---|
| 98 | betrag = 0.08 * EXP( -( radius - 2.0 * rc ) / 2.0 ) |
|---|
| 99 | ELSE |
|---|
| 100 | betrag = 0.0 |
|---|
| 101 | ENDIF |
|---|
| 102 | |
|---|
| 103 | IF ( x == 0.0 ) THEN |
|---|
| 104 | IF ( y > 0.0 ) THEN |
|---|
| 105 | alpha = pi / 2.0 |
|---|
| 106 | ELSEIF ( y < 0.0 ) THEN |
|---|
| 107 | alpha = 3.0 * pi / 2.0 |
|---|
| 108 | ENDIF |
|---|
| 109 | ELSE |
|---|
| 110 | IF ( x < 0.0 ) THEN |
|---|
| 111 | alpha = ATAN( y / x ) + pi |
|---|
| 112 | ELSE |
|---|
| 113 | IF ( y < 0.0 ) THEN |
|---|
| 114 | alpha = ATAN( y / x ) + 2.0 * pi |
|---|
| 115 | ELSE |
|---|
| 116 | alpha = ATAN( y / x ) |
|---|
| 117 | ENDIF |
|---|
| 118 | ENDIF |
|---|
| 119 | ENDIF |
|---|
| 120 | |
|---|
| 121 | vw = COS( alpha ) * betrag |
|---|
| 122 | |
|---|
| 123 | DO k = kc1, kc2 |
|---|
| 124 | v(k,j,i) = v(k,j,i) + vw |
|---|
| 125 | ENDDO |
|---|
| 126 | ENDDO |
|---|
| 127 | ENDDO |
|---|
| 128 | |
|---|
| 129 | ! |
|---|
| 130 | !-- Exchange of boundary values for the velocities. |
|---|
| 131 | CALL exchange_horiz( u, uxrp, 0 ) |
|---|
| 132 | CALL exchange_horiz( v, 0, vynp ) |
|---|
| 133 | ! |
|---|
| 134 | !-- Make velocity field nondivergent. |
|---|
| 135 | n_sor = nsor_ini |
|---|
| 136 | CALL pres |
|---|
| 137 | n_sor = nsor |
|---|
| 138 | |
|---|
| 139 | END SUBROUTINE init_rankine |
|---|
