C SUBROUTINE EMIT C C EMIT NEW VORTICES TO SATISFY BOUNDARY CONDITION. C FINISH COMPUTING PRESSURE, FORCES, ETC. C C 11/28/84 D H BAILEY MODIFIED FOR NAS KERNEL TEST C PARAMETER (NW=100, NB=5, NVM=1500) COMPLEX Z, WALL, ZCR, REFPT, EXPWKL, EXPMWK, FORCE, & UUPSTR, DUM3, EXPZ, EXPMZ, ZZ COMMON /ARRAYS/ NWALL(NB), WALL(NW,NB), RMATRX(NW*NB,NW*NB), $ ZCR(NW,NB), Z(NVM), GAMMA(NVM), REFPT(NB), RHS(NW*NB), $ FORCE(NB), RMOM(NB), CP(NW,NB), DPDS(NW,NB), EXPZ(NVM), $ EXPMZ(NVM), PSI(NW), PS(NVM) C DATA PERIOD/3./, SIG2/3./, U0/4./, MATDIM/500/, DELT/1./, $ CHORD/5./, PI/3.141592653589793/, UUPSTR/(3., 4.)/ C C STORE EXP(Z(I)) AND EXP(-Z(I)) TO REDUCE WORK IN INNER LOOP 4. C NV = 1000 PIDP = PI / PERIOD C DO 2 I = 1, NV EXPZ(I) = EXP (Z(I) * PIDP) EXPMZ(I) = 1. / EXPZ(I) 2 CONTINUE C I0 = 0 CUPST = DBLE(UUPSTR) ** 2 + IMAG(UUPSTR) ** 2 C DO 5 L = 1, NB DO 6 K = 1, NWALL(L) EXPWKL = EXP (WALL(K,L) * PIDP) EXPMWK = 1. / EXPWKL SPS = 0. DO 4 I = 1, NV DUM3 = EXPZ(I) * EXPMWK - EXPWKL * EXPMZ(I) PS(I) = GAMMA(I) * LOG (DBLE(DUM3) ** 2 + & IMAG(DUM3) ** 2 + SIG2) SPS = SPS + PS(I) 4 CONTINUE PSI(K) = IMAG(WALL(K,L) * CONJG (UUPSTR + CMPLX (0., U0))) & - SPS * 0.25 / PI 6 CONTINUE C C COMPUTE RIGHT-HAND SIDE. C DO 8 K = 1, NWALL(L) RHS(I0+K) = PSI(K) - PSI(1) 8 CONTINUE I0 = I0 + NWALL(L) 5 CONTINUE C C SOLVE SYSTEM C DO 10 I = 1, MATDIM DO 10 J = I+1, MATDIM RHS(J) = RHS(J) - RMATRX(J,I) * RHS(I) 10 CONTINUE DO 11 I = MATDIM, 1, -1 RHS(I) = RMATRX(I,I) * RHS(I) DO 11 J = 1, I-1 RHS(J) = RHS(J) - RMATRX(J,I) * RHS(I) 11 CONTINUE C C CREATE NEW VORTICES. C NOLLD = NV I0 = 0 DO 7 L = 1, NB DO 3 K = 1, NWALL(L) C C PUT THE NEW VORTEX AT THE END OF THE ARRAY. C NV = NV+1 Z(NV) = ZCR(K,L) GAMMA(NV) = RHS(I0+K) C C RECORD THE GAIN OF LINEAR AND ANGULAR MOMENTUM C FORCE(L) = FORCE(L) + GAMMA(NV) * Z(NV) RMOM(L) = RMOM(L) + GAMMA(NV) * (DBLE (Z(NV) - REFPT(L)) ** 2 & + IMAG (Z(NV) - REFPT(L)) ** 2) DPDS(K,L) = DPDS(K,L) - GAMMA(NV) 3 CONTINUE C C FILTER AND INTEGRATE PRESSURE GRADIENT TO GET PRESSURE C CP(1,L) = 0. CPM = -1E50 DO 9 K = 2, NWALL(L) CP(K,L) = CP(K-1,L) + (3. * (DPDS(K,L) + DPDS(K-1,L)) & + DPDS(1+MOD(K,NWALL(L)), L) & + DPDS(1+MOD(K+NWALL(L)-3, NWALL(L)), L)) & / (4. * DELT * CUPST) CPM = MAX (CPM, CP(K,L)) 9 CONTINUE C C NORMALIZE PRESSURE C DO 12 K = 1, NWALL(L) CP(K,L) = CP(K,L) - CPM 12 CONTINUE C C FINISH COMPUTING FORCE AND MOMENT, AS TIME RATE OF CHANGE OF LINEAR C AND ANGULAR MOMENTUM C FORCE(L) = FORCE(L) * DCMPLX (0., 2. / (DELT * CHORD * CUPST)) RMOM(L) = RMOM(L) * 2. / (DELT * CHORD ** 2 * CUPST) C I0=I0+NWALL(L) 7 CONTINUE C RETURN END