#include "dims.h" ! ! am_02/02: includine general wave (origin Kelvin wave) C SUBROUTINE DUV C **** C **** Advances u and v one time step C **** C **** where u = zonal velocity (cm/sec) C **** v = meridional velocity (cm/sec) C **** use input_module,only: iuivi,wave use init_module,only: iter use bndry_module,only: ub,ub2,uba,ubw,vb,vb2,vba,vbw, | bnd,bnd2,bnda,bndw,ci use cons_module,only: len1,len2,len3,imax,imaxh,imaxp2,kmax, | kmaxp1,expz,cor,tn,kut,gask,p0,re,dtsmooth,dtsmooth_div2, | dt,freq_semidi,shapiro,grav,dtx2inv,xmue implicit none #include "params.h" #include "fcom.h" #include "buff.h" #include "index.h" #include "phys.h" #include "vscr.h" ! real :: xmue ! COMMON/EDYVISC/XMUE(ZKMXP) C **** C **** Work space for solving the tridiagonal C **** block matrix schemes for u and v at latitude J. These C **** variables are fed to the tridiagonal block matrix C **** solver, BLKTRI. C **** ! ! Local: real :: PP(2,2,ZIMXP,ZKMXP), QQ(2,2,ZIMXP,ZKMXP), 1 RR(2,2,ZIMXP,ZKMXP), SS(2,ZIMXP,ZKMXP), BETA(2,2,ZIMXP,ZKMXP), 2 GAMMA(2,2,ZIMXP,ZKMXP),YY(2,ZIMXP,ZKMXP),XX(2,ZIMXP,ZKMXP),wper integer :: i,nuik,nvik,nflhk,nfphk,nujp2k,nujp1k,nunmk,nujm1k, | nujm2k,nvjp2k,nvjp1k,nvnmk,nvjm1k,nvjm2k,nlxxk,nlxyk,nlyxk, | nlyyk,k,ntk,nkmk,nmsk,nwk,nuk,nvk,nunpk,nvnpk,nunmnk,nvnmnk, | nxk,n C **** C **** Where: C **** PP, QQ, RR, BETA and GAMMA are (2 x 2) matrices defined C **** for k = 3/2, K-1/2, 1 and i = 3, IMAX+2, 1 C **** C **** SS, YY and XX are two component vectors similarly C **** defined. C **** C **** At level k and longitude i, the scheme is as follows: C **** C **** PP(i,k)*XX(i,k-1) + QQ(i,k)*XX(i,k) + RR(i,k)*XX(i,k+1) C **** = SS(i,k) C **** C **** BETA, GAMMA and YY are work space for intermediate C **** values used during the course of the calculation. The C **** final result is returned in XX. C **** C **** /ZTUVLB/ contains an optional lower boundary condition C **** for v as a function of latitude. C **** real :: vnlb COMMON/ZTUVLB/ VNLB(ZJMX) C **** C **** The phases of the semdiurnal, diurnal, annual tides and waves C **** are calculated as complex numbers in EXPT, EXPT2, EXPTA C **** EXPTW respectively. C **** COMPLEX EXPT,EXPT2,EXPTA,EXPTW C **** C **** Calculate u(n+1) and v(n+1) using a finite difference C **** approximation to the horizontal momentum equations. The C **** leap-frog scheme is used with vertical diffusion and C **** and advection treated implicitly to second order C **** accuracy. Since the zonal and meridional momentum C **** equations are coupled through the ion drag and Coriolis C **** terms, the system reduces to a tridiagonal block matrix C **** scheme, where (2 X 2) matrices and two component vectors C **** are used at each level. Boundary conditions for u and v C **** are required at the top and bottom of the model. C **** C **** At level k, the differential equations for u and v become: C **** C **** P(k,n)*V(k-1,n+1) + Q(k,n)*V(k,n+1) + C **** C **** R(k,n)*V(k+1,n+1) = S(k,n) (1) C **** C **** where C **** k = vertical grid index, range k = 3/2, K-1/2, 1 C **** n = time-step index C **** C **** Variables are assumed to have indices (k,n) unless C **** otherwise stated. C **** C **** P(k,n), Q(k,n), R(k,n) are (2 X 2) matrices. C **** C **** _ _ C **** | | C **** V(k,n+1) = |u(k,n+1)| (2) C **** | | C **** |v(k,n+1)| C **** |_ _| _ _ C **** | | C **** P(k,n) = -(G(k-1/2,n) + exp(-s)*w(k,n)/|1., 0.| (3) C **** (2*Ds))|0., 1.| C **** _ _ |_ _| C **** | | C **** Q(k,n) = |Q(1,1,k,n), Q(1,2,k,n)| (4) C **** | | C **** |Q(2,1,k,n), Q(2,2,k,n)| C **** |_ _| C **** C **** C **** where: C **** _ C **** Q(1,1,k,n) = G(k-1/2,n) + G(k+1/2,n) + exp(-s)* | C **** | C **** (1/(2*Dt) + lamxx(k,n)) | C **** | C **** Q(1,2,k,n) = -exp(-s)*(f+u/r*tan(theta)-lamxy) | C **** | (5) C **** Q(2,1,k,n) = exp(-s)*(f+u/r*tan(theta)-lamyx) | C **** | C **** Q(2,2,k,n) = G(k-1/2,n) + G(k+1/2,n) + exp(-s)* | C **** | C **** (1/(2*Dt) + lamyy(k,n)) _| C **** C **** _ _ C **** | | C **** R(k,n) = -(G(k+1/2,n) - exp(-s)*w(k,n)/|1., 0.| (6) C **** (2*Ds))|0., 1.| C **** |_ _| C **** C **** _ _ C **** | | C **** S(k,n) = |S(1,k,n+1)| (7) C **** | | C **** |S(2,k,n+1)| C **** |_ _| C **** C **** where: C **** _ C **** S(1,k,n) = exp(-s)*(u(n-1,k)/(2*Dt) - V.del(u) - | C **** | C **** g/(r*cos(theta))*dz/d(lamda) + lamxx*ui + | C **** | C **** lamxy*vi + FLH + GWU) | C **** (8) C **** | C **** S(2,k,n) = exp(-s)*(v(n-1,k)/(2*Dt) - V.del(v) - | C **** | C **** g/r*dz/d(theta) - lamyx*ui + lamyy*vi + FPH + GWV)| C **** _| C **** C **** G = g*KM/(P0*H*Ds**2) C **** C **** KM = kM + ke = kM + XMUE (gm/cm/sec) C **** C **** kM = vertical molecular viscosity (gm/cm/sec) C **** C **** ke = vertical eddy diffusion (gm/cm/sec) C **** C **** u = zonal neutral velocity (cm/sec) C **** C **** v = meridional neutral velocity (cm/sec) C **** C **** ui = zonal ion velocity (cm/sec) C **** C **** vi = meridional ion velocity (cm/sec) C **** C **** V.del(u) = horizontal advection of u (cm/sec**2) C **** C **** V.del(v) = horizontal advection of v (cm/sec**2) C **** C **** lamda = longitude (radians) C **** C **** theta = latitude (radians) C **** C **** r = earth's radius (cm) C **** C **** s = dimensionless vertical pressure coordinate C **** C **** Ds = grid spacing for vertical dimensionless C **** coordinate ( = 1/2) C **** C **** t = time (seconds) C **** C **** Dt = time step for leap-frog integration (seconds) C **** (typically 180 to 360 seconds) C **** _ _ C **** | | C **** |lamxx, lamxy| C **** | | = ion drag tensor (1/sec) C **** |-lamyx, lamyy| C **** |_ _| C **** C **** C **** Boundary conditions: C **** C **** Lower boundary: C **** C **** The lower boundary of the model is at (k = 1) where C **** u and v are given in one of two ways: C **** C **** a) The sum of semidiurnal, diurnal and annual tidal C **** perturbations, ub and vb C **** C **** b) Boundary conditions given by output from another C **** model e.g. CCM. The boundaries are again given C **** by ub and vb for the zonal and meridional C **** velocities respectively. C **** C **** These boundary conditions are described by: C **** C **** ub = .5*(u(1/2) + u(3/2)) C **** C **** or C **** C **** u(1/2) = 2*ub - u(3/2) C **** C **** Similarly C **** C **** v(1/2) = 2*vb - v(3/2) C **** C **** Or in vector form C **** C **** V(1/2) = 2*Vb -V(3/2) (9) C **** C **** C **** Substituting (9) in (1) for k = 3/2 gives: C **** C **** (Q(3/2,n)-P(3/2,n))*V(3/2,n+1) + R(3/2,n)*V(5/3,n+1) C **** C **** = S(3/2,n)-2.*P(3/2,n)*Vb (10) C **** C **** Upper boundary: C **** C **** The upper boundary is at level k = K, where the C **** vertical derivatives of u and v are zero. This C **** translates to: C **** C **** V(K+1/2,n+1) = V(K-1/2,n+1) (11) C **** C **** Substituting (11) in (1) for k = K-1/2 gives: C **** C **** P(K-1/2,n)*V(K-3/2,n+1) + (Q(K-1/2,n)+R(K-1/2,n))* C **** C **** V(K-1/2,n+1) = S(K-1/2,n) (12) C **** C **** Now start calculation C **** C **** Evaluate tidal boundary condition for u and v. This C **** includes semidiurnal, diurnal and annual tides. C **** C **** T1 = zonal velocity at lower boundary (cm/sec) C **** T2 = meridional velocity at lower boundary (cm/sec) C **** C **** EXPT = exp(i*nusd*t) (time varying factor for C **** semidiurnal tide) C **** C **** where i = sqrt(-1) C **** nusd = semidiurnal frequency (radians/sec) C **** t = current model time (seconds) C **** EXPT = CEXP(CI*freq_semidi*dt*ITER) C **** C **** EXPT2 = exp(i*nud*t) (time varying factor for C **** diurnal tide) C **** C **** where nud = diurnal frequency (radians/sec) C **** EXPT2 = CEXP(CI*.5*freq_semidi*dt*ITER) C **** C **** EXPTW = exp(i*nuk*t) (time varying factor for wave) C **** C **** where nuk = wave frequency (radians/sec) C **** neg. sign because Kelvin wave travels eastwards C **** wper = wave(2) ! period of wave (days) exptw = cexp(ci*.5*freq_semidi*dt/wper*iter) C **** C **** EXPTA = 1 (time varying factor for annual tide is C **** unity since it is included in the latitudinal C **** coefficients, UBA and VBA) C **** EXPTA = 1. DO I = 1,LEN1 C **** C **** Semidiurnal tide C **** T1(I) = REAL(UB(J)*BND(I)*EXPT) T2(I) = REAL(VB(J)*BND(I)*EXPT) C **** C **** Where UB = latitudinal semidiurnal tidal variation of u C **** VB = latitudinal semidiurnal tidal variation of v C **** BND = longitudinal semidiurnal tidal variation C **** of u and v C **** (UB, VB and BND are calculated in BNDRY) C **** EXPT = semidiurnal time variation C **** C **** Now add in the diurnal tide C **** T1(I) = T1(I) + REAL(UB2(J)*BND2(I)*EXPT2) T2(I) = T2(I) + REAL(VB2(J)*BND2(I)*EXPT2) C **** C **** Where UB2 = latitudinal diurnal tidal variation of u C **** VB2 = latitudinal diurnal tidal variation of v C **** BND2 = longitudinal diurnal tidal variation C **** of u and v C **** EXPT2 = diurnal time variation C **** C **** Now add in the waves C **** T1(I) = T1(I) + REAL(UBW(J)*BNDW(I)*EXPTW) T2(I) = T2(I) + REAL(VBW(J)*BNDW(I)*EXPTW) C **** C **** Where UBW = latitudinal wave variation of u C **** VBW = latitudinal wave variation of v C **** BNDW = longitudinal wave variation C **** of u and v C **** EXPTW = wave time variation C **** C **** Complete by adding annual tide. C **** T1(I) = T1(I) + REAL(UBA(J)*BNDA(I)*EXPTA) T2(I) = T2(I) + REAL(VBA(J)*BNDA(I)*EXPTA) ENDDO ! ! Report boundary conditions to stdout: ! write(6,"(/,'duv: iter=',i6,' j=',i2,' lbc for U,V')") iter,j ! write(6,"(' semidiurnal expt=',2e12.4)") expt ! write(6,"(' ub(j) =',2e12.4)") ub(j) ! write(6,"(' vb(j) =',2e12.4)") vb(j) ! write(6,"(' bnd(:) =',/,(6e12.4))") bnd ! write(6,"(' diurnal expt2=',2e12.4)") expt2 ! write(6,"(' ub2(j) =',2e12.4)") ub2(j) ! write(6,"(' vb2(j) =',2e12.4)") vb2(j) ! write(6,"(' bnd2(:)=',/,(6e12.4))") bnd2 ! write(6,"(' annual expta=',2e12.4)") expta ! write(6,"(' uba(j) =',2e12.4)") uba(j) ! write(6,"(' vba(j) =',2e12.4)") vba(j) ! write(6,"(' bnda(:)=',/,(6e12.4))") bnda ! write(6,"(' UN lbc t1=',/,(6e12.4))") t1 ! write(6,"(' VN lbc t2=',/,(6e12.4))") t2 ! C **** C **** If (IUIVI.EQ.0) set UI and VI to zero. C **** IF (IUIVI.EQ.0) THEN NUIK = NUI NVIK = NVI DO I = 1,LEN3 F(I,NUIK) = 0. F(I,NVIK) = 0. ENDDO ENDIF C **** C **** S1 = horizontal advection of u (k+1/2) C **** C **** = (u/(r*cos(theta)))*du/d(lamda) + (v/r)du/d(theta) C **** C **** (using 4th order difference approximation) C **** Note: The pointer, NU, giving the postion of the zonal C **** velocity field in the F-array, is passed to subroutine C **** ADVEC) C **** Warning: ADVEC uses work arrays S11 - S15 in /VSCR/ C **** CALL ADVEC(NU,S1) C **** C **** Similarly for v in S2 C **** CALL ADVEC(NV,S2) C **** C **** Calculate horizontal pressure forces in S14 and S15 C **** C **** S14 = g/(r*cos(theta))*dz/d(lamda) (k+1/2) C **** C **** S15 = (g/r)*dz/d(theta) (k+1/2) C **** C **** where z = geopotential height C **** C **** Warning: GLP uses S6 and S7 C **** CALL GLP(S14,S15) C **** C **** S1 = V.del(u) + g/(r*cos(theta))*dz/d(lamda) - FLH C **** C **** S2 = V.del(v) + (g/r)*dz/d(theta) - FLP C **** C **** where FLH and FPH are 4th order horizontal diffusion C **** terms. GWU and GWV are gravity wave drag terms. C **** NFLHK = NFLH NFPHK = NFPH DO I = 1,LEN2 S1(I,1) = S1(I,1) + S14(I,1) - F(I,NFLHK) S2(I,1) = S2(I,1) + S15(I,1) - F(I,NFPHK) ENDDO C **** C **** Operate Shapiro smoother on u(n-1) and v(n-1) in C **** two stages. C **** C **** First latitudinal smoothing for u with results to S6 C **** NUJP2K = NJP2 + NUNM NUJP1K = NJP1 + NUNM NUNMK = NJ + NUNM NUJM1K = NJM1 + NUNM NUJM2K = NJM2 + NUNM DO I = 1,LEN2 S6(I,1) = F(I,NUNMK) - shapiro*(F(I,NUJP2K)+F(I,NUJM2K) - 4.* 1 (F(I,NUJP1K)+F(I,NUJM1K)) + 6.*F(I,NUNMK)) ENDDO C **** C **** Now longitudinal smoothing for u with results to S7 C **** DO I = 3,LEN2-2 S7(I,1) = S6(I,1) - shapiro*(S6(I+2,1)+S6(I-2,1) - 4.* 1 (S6(I+1,1)+S6(I-1,1)) + 6.*S6(I,1)) C **** C **** S1 = V.del(u) + g/(r*cos(theta))*dz/d(lamda) - FLH - C **** 1./(2*Dt)*u(n-1) (k+1/2) C **** C **** = S1 - 1/(2*Dt)*S7 (k+1/2) C **** S1(I,1) = S1(I,1) - dtx2inv*S7(I,1) ENDDO C **** C **** Now do v C **** C **** First latitudinal smoothing for v with results to S6 C **** NVJP2K = NJP2 + NVNM NVJP1K = NJP1 + NVNM NVNMK = NJ + NVNM NVJM1K = NJM1 + NVNM NVJM2K = NJM2 + NVNM DO I = 1,LEN2 S6(I,1) = F(I,NVNMK) - shapiro*(F(I,NVJP2K)+F(I,NVJM2K) - 4.* 1 (F(I,NVJP1K)+F(I,NVJM1K)) + 6.*F(I,NVNMK)) ENDDO C **** C **** Now longitudinal smoothing for v with results to S7 C **** DO I = 3,LEN2-2 S7(I,1) = S6(I,1) - shapiro*(S6(I+2,1)+S6(I-2,1) - 4.* 1 (S6(I+1,1)+S6(I-1,1)) + 6.*S6(I,1)) C **** C **** S2 = V.del(v) + (g/r)*dz/d(theta) - FPH - C **** 1./(2*Dt)*v(n-1) (k+1/2) C **** C **** = S2 - 1/(2*Dt)*S7 (k+1/2) C **** S2(I,1) = S2(I,1) - dtx2inv*S7(I,1) ENDDO C **** C **** SS(1) = exp(-s)*(u(n-1,k)/(2*Dt) - V.del(u) - C **** C **** g/(r*cos(theta))*dz/d(lamda) + lamxx*ui + C **** C **** lamxy*vi + FLH + GWU) C **** C **** = exp(-s)*(lamxx*ui + lamxy*vi -S1) C **** C **** C **** SS(2) = exp(-s)*(v(n-1,k)/(2*Dt) - V.del(v) - C **** C **** g/r*dz/d(theta) - lamyx*ui + lamyy*vi + C **** C **** FPH + GWV) C **** C **** = exp(-s)*(lamyy*ui - lamyx*ui -S2) C **** NUIK = NUI-1 NVIK = NVI-1 NLXXK = NLXX-1 NLXYK = NLXY-1 NLYXK = NLYX-1 NLYYK = NLYY-1 DO K = 1,KMAX NUIK = NUIK+1 NVIK = NVIK+1 NLXXK = NLXXK+1 NLXYK = NLXYK+1 NLYXK = NLYXK+1 NLYYK = NLYYK+1 DO I = 1,LEN1 SS(1,I,K) = expz(K)*(.25*((F(I,NLXXK)+F(I,NLXXK+1))* 1 (F(I,NUIK)+F(I,NUIK+1))+(F(I,NLXYK)+F(I,NLXYK+1))* 2 (F(I,NVIK)+F(I,NVIK+1)))-S1(I,K)) SS(2,I,K) = expz(K)*(.25*((F(I,NLYYK)+F(I,NLYYK+1))* 1 (F(I,NVIK)+F(I,NVIK+1))-(F(I,NLYXK)+F(I,NLYXK+1))* 2 (F(I,NUIK)+F(I,NUIK+1)))-S2(I,K)) ENDDO ENDDO C **** C **** S1 = G = g*KM/(P0*H*Ds**2) C **** C **** = g**2*(kM+XMUE)*MBAR*/(P0*R*T*DS**2) C **** C **** S1 = T(k) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C **** T3 = XMUE(k) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C **** C **** Levels 2 thru KMAX C **** NTK = NJ+NT-1 DO K = 2,KMAX NTK = NTK+1 T3(K) = SQRT(XMUE(K-1)*XMUE(K)) DO I = 1,LEN1 S1(I,K) = .5*(F(I,NTK)+F(I,NTK+1)) ENDDO ENDDO C **** C **** Levels 1 and KMAXP1 C **** C **** Note: Vertical gradient of T is zero at upper C **** boundary. Value of T at lower boundary is stored in C **** F(I,NTK) at level KMAXP1. C **** NTK = NJ+NT T3(1) = SQRT(XMUE(1)**3/XMUE(2)) T3(KMAXP1) = SQRT(XMUE(KMAX)**3/XMUE(KMAX-1)) DO I = 1,LEN1 S1(I,1) = F(I,NTK+KMAX) S1(I,KMAXP1) = F(I,NTK+KMAX-1) ENDDO C **** C **** Now complete S1 = G C **** NKMK = NKM-1 NMSK = NJ+NMS-1 DO K = 1,KMAXP1 NKMK = NKMK+1 NMSK = NMSK+1 DO I = 1,LEN1 S1(I,K) = grav**2*(F(I,NKMK)+T3(K))* 1 F(I,NMSK)/(p0*gask*S1(I,K)*dz**2) ENDDO ENDDO C **** C **** S2 = exp(-s)*w(k+1/2)/(2.*Ds) C **** NWK = NJNP+NW-1 DO K = 1,KMAX NWK = NWK+1 DO I = 1,LEN1 S2(I,K) = 0.5 * expz(K) * (F(I,NWK) + F(I,NWK+1))/(2.*dz) ENDDO ENDDO C **** C **** Now evaluate: _ _ C **** | | C **** P = -(G(k-1/2,n) + exp(-s)*w(k,n)/ |1., 0.| C **** (2.*Ds)) *|0., 1.| C **** |_ _| C **** C **** and C **** _ _ C **** | | C **** R = -(G(k+1/2,n) - exp(-s)*w(k,n)/ |1., 0.| C **** (2.*Ds)) *|0., 1.| C **** |_ _| C **** for k = 3/2, K-1/2, 1 C **** DO K = 1,KMAX DO I = 1,LEN1 PP(1,1,I,K) = -(S1(I,K) + S2(I,K)) RR(1,1,I,K) = -(S1(I,K+1) - S2(I,K)) PP(2,2,I,K) = PP(1,1,I,K) RR(2,2,I,K) = RR(1,1,I,K) PP(1,2,I,K) = 0. RR(1,2,I,K) = 0. PP(2,1,I,K) = 0. RR(2,1,I,K) = 0. ENDDO ENDDO C **** C **** Now calculate Q from (5) for k = 3/2, K-1/2, 1 C **** C **** First do: C **** C **** Q(1,1,k) = G(k-1/2) + G(k+1/2) + exp(-s)*(1./(2.*Dt) + C **** lamxx(k)) C **** and C **** C **** Q(2,2,k) = G(k-1/2) + G(k+1/2) + exp(-s)*(1./(2.*Dt) + C **** lamyy(k)) C **** NLXXK = NLXX-1 NLYYK = NLYY-1 NLXYK = NLXY-1 NLYXK = NLYX-1 NUK = NJ + NU -1 DO K = 1,KMAX NLXXK = NLXXK+1 NLYYK = NLYYK+1 NLXYK = NLXYK+1 NLYXK = NLYXK+1 NUK = NUK+1 DO I = 1,LEN1 QQ(1,1,I,K) = S1(I,K) + S1(I,K+1) QQ(2,2,I,K) = QQ(1,1,I,K) QQ(1,1,I,K) = QQ(1,1,I,K) + expz(K)*(dtx2inv+.5*(F(I,NLXXK)+ 1 F(I,NLXXK+1))) QQ(2,2,I,K) = QQ(2,2,I,K) + expz(K)*(dtx2inv+.5*(F(I,NLYYK)+ 1 F(I,NLYYK+1))) C **** C **** Now do: C **** C **** Q(1,2) = -exp(-s) * (f+u/r*tan(theta)-lamxy) C **** C **** and C **** C **** Q(2,1) = exp(-s) * (f+u/r*tan(theta)-lamyx) C **** QQ(1,2,I,K) = COR(J)+F(I,NUK)/re*TN(J) QQ(2,1,I,K) = QQ(1,2,I,K) QQ(1,2,I,K) = -expz(K)*(QQ(1,2,I,K)- 1 .5*(F(I,NLXYK)+F(I,NLXYK+1))) QQ(2,1,I,K) = expz(K)*(QQ(2,1,I,K)- 1 .5*(F(I,NLYXK)+F(I,NLYXK+1))) ENDDO ENDDO C **** C **** Now set boundary conditions at k = 3/2 and k = K-1/2 C **** DO I = 1,LEN1 C **** C **** Lower boundary: C **** C **** From (10): C **** C **** Q(3/2) = Q(3/2) - P(3/2) C **** S(3/2) = S(3/2) -2.*P(3/2)*Vb C **** P(3/2) = 0. C **** QQ(1,1,I,1) = QQ(1,1,I,1) - PP(1,1,I,1) QQ(1,2,I,1) = QQ(1,2,I,1) - PP(1,2,I,1) QQ(2,1,I,1) = QQ(2,1,I,1) - PP(2,1,I,1) QQ(2,2,I,1) = QQ(2,2,I,1) - PP(2,2,I,1) SS(1,I,1) = SS(1,I,1) - 2.*(PP(1,1,I,1)*T1(I)+PP(1,2,I,1)*T2(I)) SS(2,I,1) = SS(2,I,1) - 2.*(PP(2,1,I,1)*T1(I)+PP(2,2,I,1)*T2(I)) PP(1,1,I,1) = 0. PP(1,2,I,1) = 0. PP(2,1,I,1) = 0. PP(2,2,I,1) = 0. C **** C **** Upper boundary: C **** C **** From (12): C **** C **** Q(K-1/2) = Q(K-1/2) + R(K-1/2) C **** R(K-1/2) = 0. C **** QQ(1,1,I,KMAX) = QQ(1,1,I,KMAX) + RR(1,1,I,KMAX) QQ(1,2,I,KMAX) = QQ(1,2,I,KMAX) + RR(1,2,I,KMAX) QQ(2,1,I,KMAX) = QQ(2,1,I,KMAX) + RR(2,1,I,KMAX) QQ(2,2,I,KMAX) = QQ(2,2,I,KMAX) + RR(2,2,I,KMAX) RR(1,1,I,KMAX) = 0. RR(1,2,I,KMAX) = 0. RR(2,1,I,KMAX) = 0. RR(2,2,I,KMAX) = 0. ENDDO C **** C **** Call BLKTRI to solve the tridiagonal (2 x 2) block matrix C **** for u and v. C **** CALL BLKTRI(PP,QQ,RR,SS,LEN1,3,IMAX+2,KMAXP1,1,KMAX,BETA,GAMMA, 1 YY,XX) C **** C **** Store updated u and v. C **** NUK = NJNP+NU-1 NVK = NJNP+NV-1 DO K = 1,KMAX NUK = NUK+1 NVK = NVK+1 DO I = 3,IMAX+2 F(I,NUK) = XX(1,I,K) F(I,NVK) = XX(2,I,K) ENDDO ENDDO C **** C **** Store lower boundary conditions for u and v at C **** level (KMAX+1) C **** NUK = NJNP+NU+KMAX NVK = NJNP+NV+KMAX DO I = 3,IMAX+2 F(I,NUK) =T1(I) F(I,NVK) =T2(I) ENDDO C **** C **** Fourier smoothing of u and v C **** IF(KUT(J).LT.IMAXH)THEN NUNPK = NJNP+NU NVNPK = NJNP+NV CALL FILTER(NUNPK,KMAX,KUT(J)) CALL FILTER(NVNPK,KMAX,KUT(J)) ENDIF C **** C **** Time smoothing u and v C **** NUNMK = NJ+NUNM NUK = NJ+NU NUNPK = NJNP+NU NUNMNK = NJNP+NUNM NVNMK = NJ+NVNM NVK = NJ+NV NVNPK = NJNP+NV NVNMNK = NJNP+NVNM DO I = 1,LEN2 C **** C **** u(n,smooth) = UNMN = .5*(1.-alpha)*(u(n-1)+u(n+1)) + C **** alpha*u(n) C **** F(I,NUNMNK) = dtsmooth_div2*(F(I,NUNMK)+F(I,NUNPK)) + | dtsmooth*F(I,NUK) C **** C **** and similarly for v C **** C **** F(I,NVNMNK) = dtsmooth_div2*(F(I,NVNMK)+F(I,NVNPK)) + | dtsmooth*F(I,NVK) ENDDO C **** C **** Periodic points for u and v C **** NXK = NJNP+NU-1 DO N = 1,2 DO I = 1,2 DO K = 1,2*KMAXP1 F(I,NXK+K) = F(I+IMAX,NXK+K) F(I+IMAXP2,NXK+K) = F(I+2,NXK+K) ENDDO ENDDO NXK = NJNP+NUNM-1 ENDDO RETURN END