source: trunk/SpectralModelF90/PROJECTS/trivial/rhs.f @ 6

Last change on this file since 6 was 6, checked in by xlvlod, 18 years ago

import initial from SVN_BASE_TRUNK
File size: 6.6 KB
Line 
1        subroutine rhs(F1,F2,F3,F4,F5,s2u,s2v,s2w,ds2dz,s1u,s1v,s1w,
2     *                 p_cplx,Tr1,Tr2,Tr3,Tr4,Tr5,amp,dt,myid, 
3     *                 num_dims,locnx,ny,nz,wnx,wny,wnz,
4     *                 bc_flag) 
5
6C----------------------------------------------------------------------
7C  rhs - Computes the transfer terms for the momentum and pert. density
8c        equations. N.B. all terms included except viscous/diffusive
9c        terms which are integrated analytically.
10C
11C       input: F1,F2,F3 (complex arrays of locnx by ny+1 by nz+1)
12C                        i, j, k components of the cross product of velocity 
13C                        and vorticity + rotation & buoyancy
14C              s1u, s1v, s1w   product of scalar (1) and the u, v, and w 
15C              F4        -u dot grad s1_bar
16C              F5        -u dot grad s2_bar
17c              s2u,s2v,s2w        s2*u, s2*v, s2*w
18c              ds2dz     ds2dz
19C              p_cplx    transform of d/dz s1_bar
20C              amp       amplification arrays exp(-dt*gamma*"k"^p) (see fill_wave_vectors.f)
21C              nx, ny, nz   (integers) - size of physical space data
22C              num_dims     number of dimensions 2(xz) or 3(xyz) z pos up
23C              wnx, wny, wnz (real arrays of locnx, ny+1, nz+1, respectively)
24C                            contains x, y and z wavenumbers
25C              dt         dimensionless time step
26C
27C        output: Tr1,Tr2,Tr3,Tr4 - contains the spectral representations of
28C                                  transfer terms for u, v, w and pert. density
29C              p_cplx    transform of pressure variable
30C
31C-----------------------------------------------------------------------
32
33        implicit none
34        integer locnx,ny,nz,num_dims,nyplanes,i,j,k,myid,isign
35        complex F1(locnx,ny+1,nz+1), F2(locnx,ny+1,nz+1)
36        complex F3(locnx,ny+1,nz+1), F4(locnx,ny+1,nz+1)
37        complex F5(locnx,ny+1,nz+1)
38        complex ds2dz(locnx,ny+1,nz+1)
39        complex s1u(locnx,ny+1,nz+1), s1v(locnx,ny+1,nz+1)
40        complex s1w(locnx,ny+1,nz+1),cmplxi
41        complex s2u(locnx,ny+1,nz+1), s2v(locnx,ny+1,nz+1)
42        complex s2w(locnx,ny+1,nz+1)
43        complex p_cplx(locnx,ny+1,nz+1),C
44        real amp(locnx,ny+1,nz+1,2)
45        complex Tr1(locnx,ny+1,nz+1),Tr2(locnx,ny+1,nz+1)
46        complex Tr3(locnx,ny+1,nz+1),Tr4(locnx,ny+1,nz+1)
47        complex Tr5(locnx,ny+1,nz+1)
48        real wnx(locnx),wny(ny+1),wnz(nz+1),kappa2,factor,dt
49        character*80 bc_flag
50       
51        cmplxi = (0.0,1.0)
52        if( num_dims .eq. 2 ) then
53         nyplanes=1
54        elseif( num_dims .eq. 3 ) then
55         nyplanes=ny
56        endif
57       
58        if( bc_flag .eq. 'zperiodic' ) then  ! fourier in 3d, --> d/dx --> ik
59         C= -cmplxi
60         do k = 1, nz
61          do j = 1, nyplanes
62           do i = 1,locnx
63
64            kappa2 = (wnx(i)**2) + (wny(j)**2) + (wnz(k)**2)
65            if (kappa2 .eq. 0.0) then
66             kappa2 = 1.e33
67             F1(i,j,k)=0.
68             F2(i,j,k)=0.
69             F3(i,j,k)=0.
70            endif
71            if( amp(i,j,k,2) .eq. 0 .or. wnz(k) .eq. 0 ) then
72             factor=0.0
73            else
74             factor = (-C/wnz(k))*log( amp(i,j,k,2) )/dt    ! strip -gamma*"k"^p from amp, mult by -i or 1 over m
75            endif
76
77C+
78C  scalar s1 and s2 equations: Tr4, Tr5
79C-
80            Tr4(i,j,k) = -cmplxi*wnx(i)*s1u(i,j,k)
81     *                   -cmplxi*wny(j)*s1v(i,j,k)
82     *                   -cmplxi*wnz(k)*s1w(i,j,k)
83     *                   +F4(i,j,k) + factor*p_cplx(i,j,k)   ! p_cplx has transform of d/dz s1_bar
84
85            Tr5(i,j,k) = -cmplxi*wnx(i)*s2u(i,j,k)
86     *                   -cmplxi*wny(j)*s2v(i,j,k)
87     *                   -cmplxi*wnz(k)*s2w(i,j,k)
88     *                   +F5(i,j,k) + factor*ds2dz(i,j,k)
89
90
91C+
92C  u, v & w equations: Tr1, Tr2, Tr3
93C-
94            Tr1(i,j,k) = ((wny(j)**2 + wnz(k)**2)/kappa2)*F1(i,j,k)
95     *                   - (wnx(i)*wny(j)/kappa2)*F2(i,j,k)
96     *                   - (wnx(i)*wnz(k)/kappa2)*F3(i,j,k) 
97
98            Tr2(i,j,k) = - (wnx(i)*wny(j)/kappa2)*F1(i,j,k)
99     *                   + ((wnx(i)**2 + wnz(k)**2)/kappa2)*F2(i,j,k)
100     *                   - (wny(j)*wnz(k)/kappa2)*F3(i,j,k)
101
102            Tr3(i,j,k) = - (wnx(i)*wnz(k)/kappa2)*F1(i,j,k)
103     *                   - (wny(j)*wnz(k)/kappa2)*F2(i,j,k)
104     *                   + ((wnx(i)**2 + wny(j)**2)/kappa2)*F3(i,j,k)
105
106C
107C calculate pressure 
108C
109            p_cplx(i,j,k) = -(cmplxi*wnx(i)*F1(i,j,k)+cmplxi*wny(j)*F2(i,j,k)
110     *                        +cmplxi*wnz(k)*F3(i,j,k))/kappa2                  ! overwrite w/ pressure transform
111
112           enddo
113          enddo
114         enddo
115
116        elseif( bc_flag .eq. 'zslip' ) then  ! sin/cos polarity for z derivatives
117         C= 1.0
118
119         isign=-1   ! see rhs for s1 eqn
120
121         do k = 1, nz+1
122          do j = 1, nyplanes
123           do i = 1, locnx
124
125            kappa2 = (wnx(i)**2) + (wny(j)**2) + (wnz(k)**2)
126            if (kappa2 .eq. 0.0) then
127             kappa2 = 1.e33
128             F1(i,j,k)=0.
129             F2(i,j,k)=0.
130             F3(i,j,k)=0.
131            endif
132            if( amp(i,j,k,2) .eq. 0 .or. wnz(k) .eq. 0 ) then
133             factor=0.0
134            else
135             factor = (-C/wnz(k))*log( amp(i,j,k,2) )/dt    ! strip -gamma*"k"^p from amp, mult by -i or 1 over m
136            endif
137
138C+
139C  scalar s1 and s2 equations: Tr4, Tr5
140C-
141            Tr4(i,j,k) = - cmplxi*wnx(i)*s1u(i,j,k)
142     *                   - cmplxi*wny(j)*s1v(i,j,k)
143     *                   - isign*wnz(k)*s1w(i,j,k)
144     *                   + F4(i,j,k) + factor*p_cplx(i,j,k)   ! p_cplx has transform of d/dz rho_bar on input
145
146            Tr5(i,j,k) = -cmplxi*wnx(i)*s2u(i,j,k)
147     *                   -cmplxi*wny(j)*s2v(i,j,k)
148     *                   -isign*wnz(k)*s2w(i,j,k)
149     *                   +F5(i,j,k) + factor*ds2dz(i,j,k)
150
151C+
152C  u, v & w equations: Tr1, Tr2, Tr3
153C-
154            Tr1(i,j,k) = ((wny(j)**2 + wnz(k)**2)/kappa2)*F1(i,j,k)
155     *                   - (wnx(i)*wny(j)/kappa2)*F2(i,j,k)
156     *                   + (cmplxi*wnx(i)*wnz(k)/kappa2)*F3(i,j,k) 
157
158            Tr2(i,j,k) = - (wnx(i)*wny(j)/kappa2)*F1(i,j,k)
159     *                  + ((wnx(i)**2 + wnz(k)**2)/kappa2)*F2(i,j,k)
160     *                  + (cmplxi*wny(j)*wnz(k)/kappa2)*F3(i,j,k)
161
162
163            Tr3(i,j,k) = - (cmplxi*wnx(i)*wnz(k)/kappa2)*F1(i,j,k)
164     *                   - (cmplxi*wny(j)*wnz(k)/kappa2)*F2(i,j,k)
165     *                   + ((wnx(i)**2 + wny(j)**2)/kappa2)*F3(i,j,k)
166
167C
168C calculate pressure 
169C
170            p_cplx(i,j,k) = -(cmplxi*wnx(i)*F1(i,j,k)+cmplxi*wny(j)*F2(i,j,k)
171     *                        +wnz(k)*F3(i,j,k))/kappa2                  ! overwrite w/ pressure transform
172
173           enddo
174          enddo
175         enddo
176         
177        endif
178
179        if( myid .eq. 0 ) then  ! 1,1,1 location corresponds to (k,l,m)=(0,0,0)
180         Tr1(1,1,1) = F1(1,1,1)   !  Rot*v
181         Tr2(1,1,1) = F2(1,1,1)   ! -Rot*u
182C        disallow net vertical velocity 
183         Tr3(1,1,1) = 0.
184c        set mean pressure to zero
185         p_cplx(1,1,1) = 0.
186        endif
187
188        return
189        end
Note: See TracBrowser for help on using the repository browser.