#include "dims.h" SUBROUTINE DUV use input_module,only: iuivi use bndry_module,only: ub,ub2,uba,vb,vb2,vba,bnd,bnd2,bnda,ci use cons_module,only: kut,len1,len2,len3,kmax,kmaxp1,expz,imax, | imaxp2,imaxh,cor,tn,freq_semidi,dt,freq_3m3,shapiro,dtx2inv, | p0,gask,grav,re,dtsmooth,dtsmooth_div2 implicit none C **** C **** Advances u and v one time step C **** C **** where u = zonal velocity (cm/sec) C **** v = meridional velocity (cm/sec) C **** #include "params.h" #include "fcom.h" #include "buff.h" #include "cons3m3.h" #include "index.h" #include "lowbnd.h" #include "phys.h" #include "diffk.h" #include "vscr.h" #include "bndz.h" 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),RR(2,2,ZIMXP,ZKMXP), | BETA(2,2,ZIMXP,ZKMXP),GAMMA(2,2,ZIMXP,ZKMXP), | SS(2,ZIMXP,ZKMXP), YY(2,ZIMXP,ZKMXP), XX(2,ZIMXP,ZKMXP) integer :: i,k,kk,n integer :: nuik,nvik,nujp2k,nujp1k,nunmk,nujm1k,nujm2k,nvjp2k, | nvjp1k,nvnmk,nvjm1k,nvjm2k,nlxxk,nlxyk,nlyxk,nlyyk,ntk,nmshk, | nkmk,nmsk,nwk,nuk,nunpk,nvnpk,nxk,nvk,nunmnk,nvnmnk,nfphk, | ngwvuk,ngwvvk,nflhk 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 **** /RAYFRIC/ contains values of the Rayleigh friction as a C **** function of height. C **** real :: rayk COMMON/RAYFRIC/RAYK(ZKMX) real :: raykk COMMON/GRTRAN/RAYKK(ZJMX,ZKMX) C **** C **** /EDYVISC/ contains values of the eddy viscosity as a C **** function of height. C **** C **** C **** The phases of the semdiurnal, diurnal and annual tides C **** are calculated as complex numbers in EXPT, EXPT2 and C **** EXPTA respectively. C **** C **** Complex variables for the time phase of each tidal C **** contribution C **** COMPLEX EXPT, EXPT2, EXPTA, EXPT3M3 C **** C **** Where: C **** C **** EXPT = exp(i*nu*dt*iter) C **** nu = semidiurnal frequency (radians/sec) C **** dt = timestep (seconds) C **** iter = time-step count since time zero C **** similarly, C **** C **** EXPT2 refers to the diurnal frequency C **** C **** EXPTA is the time factor for the annual tide C **** C **** EXPTM3M is the quantity needed for the (3,-3) C **** for the 2-day planetary wave C **** C **** CI = i = sqrt(-1) C **** C **** 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 + RAYK) | 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) + RAYK) _| 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 + DIFKV (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 **** C **** Evaluate tidal boundary condition for u and v. This C **** includes semidiurnal, diurnal and annual tides. Also C **** '2-day' wave. 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 **** 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 **** C EXPT2 = CEXP(CI*.5*freq_semidi*dt*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 **** C EXPTA = 1. C **** C **** EXPT3M3 = exp(i*nu3m3*t) (time varying factor for C **** (3,-3) wave) C **** where nu3m3 = frequency of '2-day' wave (radians/sec) C **** C EXPT3M3 = CEXP(CI*freq_3m3*dt*ITER) C **** C **** DO I = 1,LEN1 C **** C **** Semidiurnal tide C **** C T1(I) = REAL(UB(J)*BND(I)*EXPT) C 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 **** C T1(I) = T1(I) + REAL(UB2(J)*BND2(I)*EXPT2) C 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 **** Complete by adding annual tide. C **** C T1(I) = T1(I) + REAL(UBA(J)*BNDA(I)*EXPTA) C T2(I) = T2(I) + REAL(VBA(J)*BNDA(I)*EXPTA) C **** C **** Add in (3,-3) wave C **** C T1(I) = T1(I)+REAL(UB3M3(J)*BND3M3(I)*EXPT3M3) C T2(I) = T2(I)+REAL(VB3M3(J)*BND3M3(I)*EXPT3M3) C **** C **** where UB3M3 = latitudinal variation of (3,-3) wave C **** VB3M3 = latitudinal variation of (3,-3) wave C **** BND3M3 = longitudinal variation of (3,-3) wave C **** (TB3M3 and BND3M3 are calculated in BNDRY3M3 C **** during initialization) C **** C **** C **** Add in UBND and VBND calculated from TUVBND T1(I) = UBND(I) T2(I) = VBND(I) ENDDO ! write(6,"('duv: j=',i2,' ubnd=',/,(5e16.7))") j,ubnd ! write(6,"(12x,' vbnd=',/,(5e16.7))") vbnd 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 - GWU C **** C **** S2 = V.del(v) + (g/r)*dz/d(theta) - FLP - GWV 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 NGWVUK =NGWVU NGWVVK =NGWVV DO I = 1,LEN2 S1(I,1) = S1(I,1) + S14(I,1) - F(I,NFLHK) - F(I,NGWVUK) S2(I,1) = S2(I,1) + S15(I,1) - F(I,NFPHK) - F(I,NGWVVK) 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 **** GWU - 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 **** GWV - 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 ! 12/11/00: nuk,nvk as per kibo12: NUK = NJ + NU -1 NVK = NJ + NV -1 DO K = 1,KMAX NUIK = NUIK+1 NVIK = NVIK+1 NLXXK = NLXXK+1 NLXYK = NLXYK+1 NLYXK = NLYXK+1 NLYYK = NLYYK+1 ! 12/11/00: nuk,nvk as per kibo12: NUK = NUK+1 NVK = NVK+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)) ! 12/11/00: s8,s9 as per kibo12: S8(I,K) = (F(I,NLXXK)*(F(I,NUIK)-F(I,NUK))*1.E-2+ 1 F(I,NLXYK)*(F(I,NVIK)-F(I,NVK))*1.E-2)*86400. S9(I,K) = (F(I,NLYYK)*(F(I,NVIK)-F(I,NVK))*1.E-2- 1 F(I,NLYXK)*(F(I,NUIK)-F(I,NUK))*1.E-2)*86400. ENDDO ENDDO ! call addfsech('ALAMXX',' ',' ',F(1,NLXX),zimxp,kmaxp1,kmax,j) ! call addfsech('ALAMYY',' ',' ',F(1,NLYY),zimxp,kmaxp1,kmax,j) ! call addfsech('ALAMXY',' ',' ',F(1,NLXY),zimxp,kmaxp1,kmax,j) ! call addfsech('UII',' ',' ',F(1,NUI),zimxp,kmaxp1,kmax,j) ! call addfsech('VII',' ',' ',F(1,NVI),zimxp,kmaxp1,kmax,j) ! call addfsech('UIONM',' ',' ',S8,zimxp,kmaxp1,kmax,j) ! call addfsech('VIONM',' ',' ',S9,zimxp,kmaxp1,kmax,j) C **** C **** S1 = G = g*KM/(P0*H*Ds**2) C **** C **** = g**2*(kM+DIFKV)*MBAR*/(P0*R*T*DS**2) C **** C **** S1 = T(k) C **** S3 = rho*H**2*DIFKV (k+1/2) C **** = p0*exp(-s)*R*T*DIFKV/(MBAR*g**2) C **** NTK = NJ+NT-1 NMSHK = NMSH-1 DO K = 1,KMAX NTK = NTK+1 NMSHK = NMSHK+1 ! DO I = 1,LEN1 S3(I,K) = p0*expz(K)*gask*F(I,NTK)*DIFKV(I,K)/ 1 (F(I,NMSHK)*grav**2) ENDDO ENDDO C **** C **** Transform array S3 to full levels using logarithmic C **** interpolation. Also calculate T at full levels using C **** linear interpolation. C **** C **** Level KMAXP1 C **** C **** Note: Vertical gradient of T is zero at upper C **** boundary. C **** NTK = NJ+NT DO I = 1,LEN1 S3(I,KMAXP1) = SQRT(S3(I,KMAX)**3/S3(I,KMAX-1)) S1(I,KMAXP1) = F(I,NTK+KMAX-1) ENDDO C **** C **** Levels KMAX thru 2 C **** NTK = NJ+NT+KMAX K = KMAX+1 DO KK = 2,KMAX K = K-1 NTK = NTK-1 DO I = 1,LEN1 S3(I,K) = SQRT(S3(I,K-1)*S3(I,K)) S1(I,K) = .5*(F(I,NTK)+F(I,NTK-1)) ENDDO ENDDO C **** C **** Level 1. (Value of T at lower boundary is stored C **** in level KMAXP1) C **** NTK = NJ+NT+KMAX DO I = 1,LEN1 S3(I,1) = S3(I,1)**2/S3(I,2) S1(I,1) = F(I,NTK) 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)+S3(I,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) + RAYK(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) + RAYK(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))+RAYKK(J,K)) QQ(2,2,I,K) = QQ(2,2,I,K) + expz(K)*(dtx2inv+.5*(F(I,NLYYK)+ 1 F(I,NLYYK+1))+RAYKK(J,K)) 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,KMAXP1,KUT(J)) CALL FILTER(NVNPK,KMAXP1,KUT(J)) ENDIF C **** C **** Periodic points for u and v C **** ! 6/10/98: moved periodic points from bottom to here, so they ! can be referenced in time smoothing. 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 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,LEN3 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 RETURN END