Changeset 1001 for palm/trunk/SOURCE/init_advec.f90
 Timestamp:
 Sep 13, 2012 2:08:46 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/init_advec.f90
r667 r1001 4 4 ! Current revisions: 5 5 !  6 ! 6 ! all actions concerning upstreamsplinemethod removed 7 7 ! 8 8 ! Former revisions: … … 80 80 ENDIF 81 81 82 IF ( momentum_advec == 'upsscheme' .OR. scalar_advec == 'upsscheme' ) &83 THEN84 85 !86 ! Provide the constant parameters for the UpstreamSpline advection scheme.87 ! In x und ydirection the ShermanMorrison formula is applied88 ! (cf. Press et al, 1986 (Numerical Recipes)).89 !90 ! Allocate nonlocal arrays91 ALLOCATE( spl_z_x(0:nx), spl_z_y(0:ny), spl_tri_x(5,0:nx), &92 spl_tri_y(5,0:ny), spl_tri_zu(5,nzb:nzt+1), &93 spl_tri_zw(5,nzb:nzt+1) )94 95 !96 ! Provide diagonal elements of the tridiagonal matrices for all97 ! directions98 spl_tri_x(1,:) = 2.099 spl_tri_y(1,:) = 2.0100 spl_tri_zu(1,:) = 2.0101 spl_tri_zw(1,:) = 2.0102 103 !104 ! Elements of the cyclic tridiagonal matrix105 ! (same for all horizontal directions)106 spl_alpha = 0.5107 spl_beta = 0.5108 109 !110 ! Sub and superdiagonal elements, xdirection111 spl_tri_x(2,0:nx) = 0.5112 spl_tri_x(3,0:nx) = 0.5113 114 !115 ! mMdify the diagonal elements (ShermanMorrison)116 spl_gamma_x = spl_tri_x(1,0)117 spl_tri_x(1,0) = spl_tri_x(1,0)  spl_gamma_x118 spl_tri_x(1,nx) = spl_tri_x(1,nx)  spl_alpha * spl_beta / spl_gamma_x119 120 !121 ! Split the tridiagonal matrix for Thomas algorithm122 spl_tri_x(4,0) = spl_tri_x(1,0)123 DO i = 1, nx124 spl_tri_x(5,i) = spl_tri_x(2,i) / spl_tri_x(4,i1)125 spl_tri_x(4,i) = spl_tri_x(1,i)  spl_tri_x(5,i) * spl_tri_x(3,i1)126 ENDDO127 128 !129 ! Allocate arrays required locally130 ALLOCATE( temp(0:nx), spl_u(0:nx) )131 132 !133 ! Provide "corrective vector", xdirection134 spl_u(0) = spl_gamma_x135 spl_u(1:nx1) = 0.0136 spl_u(nx) = spl_alpha137 138 !139 ! Solve the system of equations for the corrective vector140 ! (ShermanMorrison)141 temp(0) = spl_u(0)142 DO i = 1, nx143 temp(i) = spl_u(i)  spl_tri_x(5,i) * temp(i1)144 ENDDO145 spl_z_x(nx) = temp(nx) / spl_tri_x(4,nx)146 DO i = nx1, 0, 1147 spl_z_x(i) = ( temp(i)  spl_tri_x(3,i) * spl_z_x(i+1) ) / &148 spl_tri_x(4,i)149 ENDDO150 151 !152 ! Deallocate local arrays, for they are allocated in a different way for153 ! operations in ydirection154 DEALLOCATE( temp, spl_u )155 156 !157 ! Provide sub and superdiagonal elements, ydirection158 spl_tri_y(2,0:ny) = 0.5159 spl_tri_y(3,0:ny) = 0.5160 161 !162 ! Modify the diagonal elements (ShermanMorrison)163 spl_gamma_y = spl_tri_y(1,0)164 spl_tri_y(1,0) = spl_tri_y(1,0)  spl_gamma_y165 spl_tri_y(1,ny) = spl_tri_y(1,ny)  spl_alpha * spl_beta / spl_gamma_y166 167 !168 ! Split the tridiagonal matrix for Thomas algorithm169 spl_tri_y(4,0) = spl_tri_y(1,0)170 DO j = 1, ny171 spl_tri_y(5,j) = spl_tri_y(2,j) / spl_tri_y(4,j1)172 spl_tri_y(4,j) = spl_tri_y(1,j)  spl_tri_y(5,j) * spl_tri_y(3,j1)173 ENDDO174 175 !176 ! Allocate arrays required locally177 ALLOCATE( temp(0:ny), spl_u(0:ny) )178 179 !180 ! Provide "corrective vector", ydirection181 spl_u(0) = spl_gamma_y182 spl_u(1:ny1) = 0.0183 spl_u(ny) = spl_alpha184 185 !186 ! Solve the system of equations for the corrective vector187 ! (ShermanMorrison)188 temp = 0.0189 spl_z_y = 0.0190 temp(0) = spl_u(0)191 DO j = 1, ny192 temp(j) = spl_u(j)  spl_tri_y(5,j) * temp(j1)193 ENDDO194 spl_z_y(ny) = temp(ny) / spl_tri_y(4,ny)195 DO j = ny1, 0, 1196 spl_z_y(j) = ( temp(j)  spl_tri_y(3,j) * spl_z_y(j+1) ) / &197 spl_tri_y(4,j)198 ENDDO199 200 !201 ! deallocate local arrays, for they are no longer required202 DEALLOCATE( temp, spl_u )203 204 !205 ! provide sub and superdiagonal elements, zdirection206 spl_tri_zu(2,nzb) = 0.0207 spl_tri_zu(2,nzt+1) = 1.0208 spl_tri_zw(2,nzb) = 0.0209 spl_tri_zw(2,nzt+1) = 1.0210 211 spl_tri_zu(3,nzb) = 1.0212 spl_tri_zu(3,nzt+1) = 0.0213 spl_tri_zw(3,nzb) = 1.0214 spl_tri_zw(3,nzt+1) = 0.0215 216 DO k = nzb+1, nzt217 spl_tri_zu(2,k) = dzu(k) / ( dzu(k) + dzu(k+1) )218 spl_tri_zw(2,k) = dzw(k) / ( dzw(k) + dzw(k+1) )219 spl_tri_zu(3,k) = 1.0  spl_tri_zu(2,k)220 spl_tri_zw(3,k) = 1.0  spl_tri_zw(2,k)221 ENDDO222 223 spl_tri_zu(4,nzb) = spl_tri_zu(1,nzb)224 spl_tri_zw(4,nzb) = spl_tri_zw(1,nzb)225 DO k = nzb+1, nzt+1226 spl_tri_zu(5,k) = spl_tri_zu(2,k) / spl_tri_zu(4,k1)227 spl_tri_zw(5,k) = spl_tri_zw(2,k) / spl_tri_zw(4,k1)228 spl_tri_zu(4,k) = spl_tri_zu(1,k)  spl_tri_zu(5,k) * &229 spl_tri_zu(3,k1)230 spl_tri_zw(4,k) = spl_tri_zw(1,k)  spl_tri_zw(5,k) * &231 spl_tri_zw(3,k1)232 ENDDO233 234 ENDIF235 236 82 END SUBROUTINE init_advec
Note: See TracChangeset
for help on using the changeset viewer.