SUBROUTINE SICDKD (D,T,Q,X,Y,Z,AM,BM,CM,TOV,SV,DT,ZE,W,GZ,DISPK) C D, T, Q, ARE SPECTRAL COMPONENTS AT (T-1) C X, Y, Z ARE SPECTRAL COMPONENTS OF NONLIN. TERMS DIMENSION D( 289 , 7 ), T( 289 , 7 ), 1 Q( 289 ), X( 289 , 7 ), 2 Y( 289 , 7 ), Z( 289 ), AM( 7 , 7 ), BM( 7 , 7 ), 3 CM( 7 , 7 ), TOV( 7 ), SV( 7 ) DIMENSION ZE( 289 , 7 ), W( 289 , 7 ) ,GZ( 289 ) COMPLEX ZE, W C VARIABLES USED ONLY IN SUB. DIMENSION RIMDIS( 7 , 7 ) DIMENSION RIM( 7 , 7 ), DM( 7 , 7 ), 1 EM( 7 , 7 ), LV( 7 ), MV( 7 ), XV( 7 ), 2 YV( 7 ), SAV( 7 ), SBV( 7 ), 3 RHSV( 7 ), DV( 7 ), TV( 7 ) DIMENSION TOVRA( 7 ) COMPLEX T, D, Q, X, Y, Z,GZ COMPLEX XV, YV, SAV, SBV, RHSV, DV, TV, SUM JCAP = 16 LEVS = 7 JCAP1 = JCAP + 1 C DT2 = 2. * DT DT2 = DT + DT C TOV = R*TO/Q**2 C SET UNIT MATRIX RIM DO 1 J=1,LEVS C TOVRA(J)= 287.05 *TOV(J)/( 6370000. * 6370000. ) TOVRA(J)=( 287.05 /( 6370000. * 6370000. ))*TOV(J) DO 1 I=1,LEVS RIMDIS(I,J)=0. RIM(I,J) = 0. 1 CONTINUE DO 2 I=1,LEVS RIM(I,I) = 1. RIMDIS(I,I)=DT2*DISPK +1. 2 CONTINUE CCC 1000 IS (L,N) SPECTRAL COMPONENT LOOP LL1=1 JMN=2* 16 +1 DO 1000 NN=1,JMN N=NN-1 C COMPUTE MATRICES DM, EM FROM CM DO 4 J=1,LEVS DO 4 I=1,LEVS TEMP = N * (N+1.) * CM(I,J) C DM(I,J) = RIM(I,J) - TEMP DM(I,J)=RIMDIS(I,J)-TEMP EM(I,J) = RIM(I,J) + TEMP 4 CONTINUE CALL IMINV (DM, LEVS, DET, LV, MV) C DM NOW CONTAINS ITS INVERSE. LL2=NN IF(N- 16 )900,900,800 800 LL1=LL1+1 LL2= 17 900 CONTINUE IUP=0 DO 1000 LL=LL1,LL2 L=LL2-LL IUP=IUP+1 IF(N.GT. 16 )L=LL2-IUP IN =L* 17 +LL C SET COLUMN VECTORS XV, YV, DV, TV DO 5 LEV=1,LEVS DV(LEV) = D(IN,LEV) TV(LEV) = T(IN,LEV) XV(LEV) = X(IN,LEV) YV(LEV) = Y(IN,LEV) 5 CONTINUE C BEGIN CONSTRUCTING RHSV (FROM R TO L) DO 6 LEV=1,LEVS C RHSV(LEV)=(DT*Z(IN)+Q(IN))*TOVRA(LEV) RHSV(LEV) = CMPLX((REAL(Z(IN))*DT + REAL(Q(IN)))*TOVRA(LEV), 1 (AIMAG(Z(IN))*DT + AIMAG(Q(IN)))*TOVRA(LEV)) 6 CONTINUE DO 7 LEV=1,LEVS C SAV(LEV) = DT * YV(LEV) + TV(LEV) SAV(LEV) = CMPLX(DT*REAL(YV(LEV)) + REAL(TV(LEV)), 1 DT*AIMAG(YV(LEV)) + AIMAG(TV(LEV))) 7 CONTINUE C MULTIPLY MATRIX AM BY SAV PUT IN SBV CALL MATVEC (AM, SAV, SBV) DO 10 LEV=1,LEVS C RHSV(LEV) = N * (N+1.) * (RHSV(LEV) + SBV(LEV)) RHSV(LEV) = CMPLX((REAL(RHSV(LEV)) + REAL(SBV(LEV)))*(N*(N+1)), 1 (AIMAG(RHSV(LEV)) + AIMAG(SBV(LEV)))*(N*(N+1.))) 10 CONTINUE DO 11 LEV=1,LEVS C RHSV(LEV) = DT2 * (XV(LEV) + RHSV(LEV)) RHSV(LEV) = CMPLX((REAL(XV(LEV)) + REAL(RHSV(LEV)))*DT2, 1 (AIMAG(XV(LEV)) + AIMAG(RHSV(LEV)))*DT2) 11 CONTINUE C MULTIPLY EM MATRIX BY DV(T-1) PUT IN SAV CALL MATVEC (EM, DV, SAV) DO 14 LEV=1,LEVS RHSV(LEV) = RHSV(LEV) + SAV(LEV) 14 CONTINUE C MULTIPLY RHSV BY DM INVERSE (NOW IN DM) CALL MATVEC (DM, RHSV, SAV) C SAV IS NOW D VECTOR AT T+1 C STORE D AT T+1 IN X DO 15 LEV=1,LEVS X(IN,LEV) = SAV(LEV) 15 CONTINUE C STEP T, Q IN TIME. C SET SAV = D(T+1) + D(T-1) DO 16 LEV=1,LEVS SAV(LEV) = SAV(LEV) + DV(LEV) 16 CONTINUE C MULTIPLY BM BY SAV PUT IN SBV CALL MATVEC (BM, SAV, SBV) DO 20 LEV=1,LEVS C Y(IN,LEV) = TV(LEV) + DT * (2. * YV(LEV) + SBV(LEV)) Y(IN,LEV) = CMPLX((REAL(YV(LEV))+REAL(YV(LEV))+REAL(SBV(LEV)))*DT+ 1 REAL(TV(LEV)), 2 (AIMAG(YV(LEV))+AIMAG(YV(LEV))+AIMAG(SBV(LEV)))* 3 DT + AIMAG(TV(LEV))) 20 CONTINUE C MULTIPLY S TRANSPOSE BY SAV SUM = CMPLX(0.,0.) DO 21 LEV=1,LEVS C SUM = SUM + SV(LEV) * SAV(LEV) SUM = CMPLX(REAL(SUM) + SV(LEV)*REAL(SAV(LEV)), 1 AIMAG(SUM) + SV(LEV)*AIMAG(SAV(LEV))) 21 CONTINUE C Z(IN) = Q(IN) + DT * (2. * Z(IN) + SUM) Z(IN) = CMPLX(REAL(Q(IN))+(REAL(Z(IN))+REAL(Z(IN))+REAL(SUM))*DT, 1 AIMAG(Q(IN))+(AIMAG(Z(IN))+AIMAG(Z(IN))+AIMAG(SUM))*DT) 1000 CONTINUE C COMPUTE VORTICITY EXPLICITLY DO 30 LEV=1, 7 DO 30 I=1, 289 C W(I,LEV) = ZE(I,LEV) + 2.*DT*W(I,LEV) W(I,LEV) = CMPLX(REAL(ZE(I,LEV)) + REAL(W(I,LEV))*DT2, 1 AIMAG(ZE(I,LEV)) + AIMAG(W(I,LEV))*DT2) 30 CONTINUE C(DEBUGLINE) VERSION CREATED JULY 8 1982 BY J.SELA,NMC. RETURN END