[1] | 1 | SUBROUTINE init_advec |
---|
| 2 | |
---|
| 3 | !------------------------------------------------------------------------------! |
---|
| 4 | ! Actual revisions: |
---|
| 5 | ! ----------------- |
---|
| 6 | ! |
---|
| 7 | ! |
---|
| 8 | ! Former revisions: |
---|
| 9 | ! ----------------- |
---|
| 10 | ! $Log: init_advec.f90,v $ |
---|
| 11 | ! Revision 1.6 2004/04/30 11:59:31 raasch |
---|
| 12 | ! impulse_advec renamed momentum_advec |
---|
| 13 | ! |
---|
| 14 | ! Revision 1.5 2001/03/30 07:25:53 raasch |
---|
| 15 | ! Translation of remaining German identifiers (variables, subroutines, etc.) |
---|
| 16 | ! |
---|
| 17 | ! Revision 1.4 2001/01/22 06:57:43 raasch |
---|
| 18 | ! Module test_variables removed |
---|
| 19 | ! |
---|
| 20 | ! Revision 1.3 2000/12/21 16:12:52 letzel |
---|
| 21 | ! All comments translated into English. |
---|
| 22 | ! |
---|
| 23 | ! Revision 1.2 1999/02/17 09:18:45 raasch |
---|
| 24 | ! Fehlerkorrektur, Verfahren muessen unabhaengig voneinander initialisierbar |
---|
| 25 | ! sein |
---|
| 26 | ! |
---|
| 27 | ! Revision 1.1 1999/02/05 09:07:38 raasch |
---|
| 28 | ! Initial revision |
---|
| 29 | ! |
---|
| 30 | ! |
---|
| 31 | ! Description: |
---|
| 32 | ! ------------ |
---|
| 33 | ! Initialize constant coefficients and parameters for certain advection schemes. |
---|
| 34 | !------------------------------------------------------------------------------! |
---|
| 35 | |
---|
| 36 | USE advection |
---|
| 37 | USE arrays_3d |
---|
| 38 | USE indices |
---|
| 39 | USE control_parameters |
---|
| 40 | |
---|
| 41 | IMPLICIT NONE |
---|
| 42 | |
---|
| 43 | INTEGER :: i, intervals, j, k |
---|
| 44 | REAL :: delt, dn, dnneu, ex1, ex2, ex3, ex4, ex5, ex6, spl_alpha, & |
---|
| 45 | spl_beta, sterm |
---|
| 46 | REAL, DIMENSION(:), ALLOCATABLE :: spl_u, temp |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | IF ( scalar_advec == 'bc-scheme' ) THEN |
---|
| 50 | |
---|
| 51 | ! |
---|
| 52 | !-- Compute exponential coefficients for the Bott-Chlond scheme |
---|
| 53 | intervals = 1000 |
---|
| 54 | ALLOCATE( aex(intervals), bex(intervals), dex(intervals), eex(intervals) ) |
---|
| 55 | |
---|
| 56 | delt = 1.0 / REAL( intervals ) |
---|
| 57 | sterm = delt * 0.5 |
---|
| 58 | |
---|
| 59 | DO i = 1, intervals |
---|
| 60 | |
---|
| 61 | IF ( sterm > 0.5 ) THEN |
---|
| 62 | dn = -5.0 |
---|
| 63 | ELSE |
---|
| 64 | dn = 5.0 |
---|
| 65 | ENDIF |
---|
| 66 | |
---|
| 67 | DO j = 1, 15 |
---|
| 68 | ex1 = dn * EXP( -dn ) - EXP( 0.5 * dn ) + EXP( -0.5 * dn ) |
---|
| 69 | ex2 = EXP( dn ) - EXP( -dn ) |
---|
| 70 | ex3 = EXP( -dn ) * ( 1.0 - dn ) - 0.5 * EXP( 0.5 * dn ) & |
---|
| 71 | - 0.5 * EXP( -0.5 * dn ) |
---|
| 72 | ex4 = EXP( dn ) + EXP( -dn ) |
---|
| 73 | ex5 = dn * sterm + ex1 / ex2 |
---|
| 74 | ex6 = sterm + ( ex3 * ex2 - ex4 * ex1 ) / ( ex2 * ex2 ) |
---|
| 75 | dnneu = dn - ex5 / ex6 |
---|
| 76 | dn = dnneu |
---|
| 77 | ENDDO |
---|
| 78 | |
---|
| 79 | IF ( sterm < 0.5 ) dn = MAX( 2.95E-2, dn ) |
---|
| 80 | IF ( sterm > 0.5 ) dn = MIN( -2.95E-2, dn ) |
---|
| 81 | ex1 = EXP( -dn ) |
---|
| 82 | ex2 = EXP( dn ) - ex1 |
---|
| 83 | aex(i) = -ex1 / ex2 |
---|
| 84 | bex(i) = 1.0 / ex2 |
---|
| 85 | dex(i) = dn |
---|
| 86 | eex(i) = EXP( dex(i) * 0.5 ) |
---|
| 87 | sterm = sterm + delt |
---|
| 88 | |
---|
| 89 | ENDDO |
---|
| 90 | |
---|
| 91 | ENDIF |
---|
| 92 | |
---|
| 93 | IF ( momentum_advec == 'ups-scheme' .OR. scalar_advec == 'ups-scheme' ) & |
---|
| 94 | THEN |
---|
| 95 | |
---|
| 96 | ! |
---|
| 97 | !-- Provide the constant parameters for the Upstream-Spline advection scheme. |
---|
| 98 | !-- In x- und y-direction the Sherman-Morrison formula is applied |
---|
| 99 | !-- (cf. Press et al, 1986 (Numerical Recipes)). |
---|
| 100 | ! |
---|
| 101 | !-- Allocate nonlocal arrays |
---|
| 102 | ALLOCATE( spl_z_x(0:nx), spl_z_y(0:ny), spl_tri_x(5,0:nx), & |
---|
| 103 | spl_tri_y(5,0:ny), spl_tri_zu(5,nzb:nzt+1), & |
---|
| 104 | spl_tri_zw(5,nzb:nzt+1) ) |
---|
| 105 | |
---|
| 106 | ! |
---|
| 107 | !-- Provide diagonal elements of the tridiagonal matrices for all |
---|
| 108 | !-- directions |
---|
| 109 | spl_tri_x(1,:) = 2.0 |
---|
| 110 | spl_tri_y(1,:) = 2.0 |
---|
| 111 | spl_tri_zu(1,:) = 2.0 |
---|
| 112 | spl_tri_zw(1,:) = 2.0 |
---|
| 113 | |
---|
| 114 | ! |
---|
| 115 | !-- Elements of the cyclic tridiagonal matrix |
---|
| 116 | !-- (same for all horizontal directions) |
---|
| 117 | spl_alpha = 0.5 |
---|
| 118 | spl_beta = 0.5 |
---|
| 119 | |
---|
| 120 | ! |
---|
| 121 | !-- Sub- and superdiagonal elements, x-direction |
---|
| 122 | spl_tri_x(2,0:nx) = 0.5 |
---|
| 123 | spl_tri_x(3,0:nx) = 0.5 |
---|
| 124 | |
---|
| 125 | ! |
---|
| 126 | !-- mMdify the diagonal elements (Sherman-Morrison) |
---|
| 127 | spl_gamma_x = -spl_tri_x(1,0) |
---|
| 128 | spl_tri_x(1,0) = spl_tri_x(1,0) - spl_gamma_x |
---|
| 129 | spl_tri_x(1,nx) = spl_tri_x(1,nx) - spl_alpha * spl_beta / spl_gamma_x |
---|
| 130 | |
---|
| 131 | ! |
---|
| 132 | !-- Split the tridiagonal matrix for Thomas algorithm |
---|
| 133 | spl_tri_x(4,0) = spl_tri_x(1,0) |
---|
| 134 | DO i = 1, nx |
---|
| 135 | spl_tri_x(5,i) = spl_tri_x(2,i) / spl_tri_x(4,i-1) |
---|
| 136 | spl_tri_x(4,i) = spl_tri_x(1,i) - spl_tri_x(5,i) * spl_tri_x(3,i-1) |
---|
| 137 | ENDDO |
---|
| 138 | |
---|
| 139 | ! |
---|
| 140 | !-- Allocate arrays required locally |
---|
| 141 | ALLOCATE( temp(0:nx), spl_u(0:nx) ) |
---|
| 142 | |
---|
| 143 | ! |
---|
| 144 | !-- Provide "corrective vector", x-direction |
---|
| 145 | spl_u(0) = spl_gamma_x |
---|
| 146 | spl_u(1:nx-1) = 0.0 |
---|
| 147 | spl_u(nx) = spl_alpha |
---|
| 148 | |
---|
| 149 | ! |
---|
| 150 | !-- Solve the system of equations for the corrective vector |
---|
| 151 | !-- (Sherman-Morrison) |
---|
| 152 | temp(0) = spl_u(0) |
---|
| 153 | DO i = 1, nx |
---|
| 154 | temp(i) = spl_u(i) - spl_tri_x(5,i) * temp(i-1) |
---|
| 155 | ENDDO |
---|
| 156 | spl_z_x(nx) = temp(nx) / spl_tri_x(4,nx) |
---|
| 157 | DO i = nx-1, 0, -1 |
---|
| 158 | spl_z_x(i) = ( temp(i) - spl_tri_x(3,i) * spl_z_x(i+1) ) / & |
---|
| 159 | spl_tri_x(4,i) |
---|
| 160 | ENDDO |
---|
| 161 | |
---|
| 162 | ! |
---|
| 163 | !-- Deallocate local arrays, for they are allocated in a different way for |
---|
| 164 | !-- operations in y-direction |
---|
| 165 | DEALLOCATE( temp, spl_u ) |
---|
| 166 | |
---|
| 167 | ! |
---|
| 168 | !-- Provide sub- and superdiagonal elements, y-direction |
---|
| 169 | spl_tri_y(2,0:ny) = 0.5 |
---|
| 170 | spl_tri_y(3,0:ny) = 0.5 |
---|
| 171 | |
---|
| 172 | ! |
---|
| 173 | !-- Modify the diagonal elements (Sherman-Morrison) |
---|
| 174 | spl_gamma_y = -spl_tri_y(1,0) |
---|
| 175 | spl_tri_y(1,0) = spl_tri_y(1,0) - spl_gamma_y |
---|
| 176 | spl_tri_y(1,ny) = spl_tri_y(1,ny) - spl_alpha * spl_beta / spl_gamma_y |
---|
| 177 | |
---|
| 178 | ! |
---|
| 179 | !-- Split the tridiagonal matrix for Thomas algorithm |
---|
| 180 | spl_tri_y(4,0) = spl_tri_y(1,0) |
---|
| 181 | DO j = 1, ny |
---|
| 182 | spl_tri_y(5,j) = spl_tri_y(2,j) / spl_tri_y(4,j-1) |
---|
| 183 | spl_tri_y(4,j) = spl_tri_y(1,j) - spl_tri_y(5,j) * spl_tri_y(3,j-1) |
---|
| 184 | ENDDO |
---|
| 185 | |
---|
| 186 | ! |
---|
| 187 | !-- Allocate arrays required locally |
---|
| 188 | ALLOCATE( temp(0:ny), spl_u(0:ny) ) |
---|
| 189 | |
---|
| 190 | ! |
---|
| 191 | !-- Provide "corrective vector", y-direction |
---|
| 192 | spl_u(0) = spl_gamma_y |
---|
| 193 | spl_u(1:ny-1) = 0.0 |
---|
| 194 | spl_u(ny) = spl_alpha |
---|
| 195 | |
---|
| 196 | ! |
---|
| 197 | !-- Solve the system of equations for the corrective vector |
---|
| 198 | !-- (Sherman-Morrison) |
---|
| 199 | temp = 0.0 |
---|
| 200 | spl_z_y = 0.0 |
---|
| 201 | temp(0) = spl_u(0) |
---|
| 202 | DO j = 1, ny |
---|
| 203 | temp(j) = spl_u(j) - spl_tri_y(5,j) * temp(j-1) |
---|
| 204 | ENDDO |
---|
| 205 | spl_z_y(ny) = temp(ny) / spl_tri_y(4,ny) |
---|
| 206 | DO j = ny-1, 0, -1 |
---|
| 207 | spl_z_y(j) = ( temp(j) - spl_tri_y(3,j) * spl_z_y(j+1) ) / & |
---|
| 208 | spl_tri_y(4,j) |
---|
| 209 | ENDDO |
---|
| 210 | |
---|
| 211 | ! |
---|
| 212 | !-- deallocate local arrays, for they are no longer required |
---|
| 213 | DEALLOCATE( temp, spl_u ) |
---|
| 214 | |
---|
| 215 | ! |
---|
| 216 | !-- provide sub- and superdiagonal elements, z-direction |
---|
| 217 | spl_tri_zu(2,nzb) = 0.0 |
---|
| 218 | spl_tri_zu(2,nzt+1) = 1.0 |
---|
| 219 | spl_tri_zw(2,nzb) = 0.0 |
---|
| 220 | spl_tri_zw(2,nzt+1) = 1.0 |
---|
| 221 | |
---|
| 222 | spl_tri_zu(3,nzb) = 1.0 |
---|
| 223 | spl_tri_zu(3,nzt+1) = 0.0 |
---|
| 224 | spl_tri_zw(3,nzb) = 1.0 |
---|
| 225 | spl_tri_zw(3,nzt+1) = 0.0 |
---|
| 226 | |
---|
| 227 | DO k = nzb+1, nzt |
---|
| 228 | spl_tri_zu(2,k) = dzu(k) / ( dzu(k) + dzu(k+1) ) |
---|
| 229 | spl_tri_zw(2,k) = dzw(k) / ( dzw(k) + dzw(k+1) ) |
---|
| 230 | spl_tri_zu(3,k) = 1.0 - spl_tri_zu(2,k) |
---|
| 231 | spl_tri_zw(3,k) = 1.0 - spl_tri_zw(2,k) |
---|
| 232 | ENDDO |
---|
| 233 | |
---|
| 234 | spl_tri_zu(4,nzb) = spl_tri_zu(1,nzb) |
---|
| 235 | spl_tri_zw(4,nzb) = spl_tri_zw(1,nzb) |
---|
| 236 | DO k = nzb+1, nzt+1 |
---|
| 237 | spl_tri_zu(5,k) = spl_tri_zu(2,k) / spl_tri_zu(4,k-1) |
---|
| 238 | spl_tri_zw(5,k) = spl_tri_zw(2,k) / spl_tri_zw(4,k-1) |
---|
| 239 | spl_tri_zu(4,k) = spl_tri_zu(1,k) - spl_tri_zu(5,k) * & |
---|
| 240 | spl_tri_zu(3,k-1) |
---|
| 241 | spl_tri_zw(4,k) = spl_tri_zw(1,k) - spl_tri_zw(5,k) * & |
---|
| 242 | spl_tri_zw(3,k-1) |
---|
| 243 | ENDDO |
---|
| 244 | |
---|
| 245 | ENDIF |
---|
| 246 | |
---|
| 247 | END SUBROUTINE init_advec |
---|