m0n0lithic

SUBROUTINE VORTX1

Jun 4th, 2014
238
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.       SUBROUTINE VORTX1
  2.  
  3. c     Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000
  4.  
  5.       INCLUDE 'COMML.FOR'
  6.       epsp=0.0001
  7.       MOL=0
  8.       DO 19 I=1,M1
  9.       DO 19 J=1,N1
  10.       IF(KT.EQ.1)PSI(J,I)=PSI0(J,I)
  11.       IF(H(J,I).EQ.PF)GOTO 19
  12.       MOL=MOL+1
  13.    19 CONTINUE
  14.       TETA=1.
  15.       L=0  
  16.    35 L=L+1
  17.       IF(L.GE.2000000)GOTO 325
  18.  
  19.       CALL BOUND
  20.  
  21.       DO 778 I=2,M
  22.       DO 778 J=2,N
  23.       IF(H(J,I).EQ.PF)GOTO 778
  24. C     IF(H(J+1,I).EQ.PF.OR.H(J,I+1).EQ.PF.OR.H(J,I-1).EQ.PF
  25. C    &.OR.H(J-1,I).EQ.PF)GOTO 778
  26. C     IF(H(J+1,I+1).EQ.PF.OR.H(J-1,I+1).EQ.PF.OR.H(J+1,I-1).EQ.PF
  27. C    &.OR.H(J-1,I-1).EQ.PF)GOTO 778
  28.       DDT=(AFV(J,I)+AFV(J-1,I)-AFV(J,I-1)-AFV(J-1,I-1)
  29.      &-AFU(J,I)-AFU(J,I-1)+AFU(J-1,I)+AFU(J-1,I-1))
  30.  
  31.       A1=RK(J,I)*HL1(J,I)/HL2(J,I)/H1(J,I)
  32.       A4=RK(J,I-1)*HL1(J,I-1)/HL2(J,I-1)/H1(J,I-1)
  33.       A8=RK(J-1,I)*HL1(J-1,I)/HL2(J-1,I)/H1(J-1,I)
  34.       A9=RK(J-1,I-1)*HL1(J-1,I-1)/HL2(J-1,I-1)/H1(J-1,I-1)
  35.       B4=YL*HH1(J+1,I)/HH2(J+1,I)/H(J+1,I)
  36.       B8=YL*HH1(J-1,I)/HH2(J-1,I)/H(J-1,I)
  37.       G1=RK(J,I)*HL2(J,I)/HL1(J,I)/H1(J,I)
  38.       G2=RK(J-1,I)*HL2(J-1,I)/HL1(J-1,I)/H1(J-1,I)
  39.       G6=RK(J,I-1)*HL2(J,I-1)/HL1(J,I-1)/H1(J,I-1)
  40.       G7=RK(J-1,I-1)*HL2(J-1,I-1)/HL1(J-1,I-1)/H1(J-1,I-1)
  41.       O2=YL*HH2(J,I+1)/HH1(J,I+1)/H(J,I+1)
  42.       O6=YL*HH2(J,I-1)/HH1(J,I-1)/H(J,I-1)
  43.  
  44.       PSI1=((G2+G1)*PSI0(J,I+1)+(G7+G6)*PSI(J,I-1)+
  45.      & (A4+A1)*PSI0(J+1,I)+(A9+A8)*PSI(J-1,I)+
  46.      & 0.5*(B4*(PSI0(J+1,I+1)-PSI0(J+1,I-1))-B8*(PSI0(J-1,I+1)-
  47.      & PSI(J-1,I-1))-O2*(PSI0(J+1,I+1)-PSI0(J-1,I+1))+
  48.      & O6*(PSI0(J+1,I-1)-PSI(J-1,I-1)))-HAG*DDT)/         ! DOESN'T DIVIDE BY 2 AND ONE HAG!
  49.      & (G1+G2+G6+G7+A1+A4+A8+A9)
  50.  
  51.       PSI(J,I)=TETA*PSI1+(1.-TETA)*PSI0(J,I)
  52.  778  CONTINUE
  53.  
  54.       CALL VELCTY
  55.       CALL FRIC
  56.       CALL RHS
  57.  
  58.       if(L.eq.1)go to 35
  59.       KAP=0
  60.       DO 97 I=1,M1
  61.       DO 97 J=1,N1
  62.       IF(H(J,I).EQ.PF)GOTO 97
  63.       IF(PSI0(J,I).EQ.0..AND.PSI(J,I).EQ.0.)GOTO 322
  64.       IF(ABS((PSI(J,I)-PSI0(J,I))/(PSI0(J,I)+0.1E-09)).GT.EPSP)GOTO 97
  65.  322  KAP=KAP+1
  66.   97  CONTINUE
  67.       DO 2 I=1,M1
  68.       DO 2 J=1,N1
  69.       IF(H(J,I).EQ.PF)GOTO 2
  70.       PSI0(J,I)=PSI(J,I)
  71.    2  CONTINUE
  72.       IF(KAP.EQ.MOL)GOTO 444
  73.       IF(L/ITER1*ITER1.EQ.L)WRITE(*,*)
  74.      &'Acercamiento esta en ',KAP/FLOAT(MOL)*100,' %'
  75.       IF(KAP/FLOAT(MOL)*100..GE.PERC)GO TO 444
  76.       GOTO 35
  77.  444  CONTINUE
  78. CCC   write(*,*)l,kap,mol
  79.       RETURN
  80.  325  WRITE(*,764)
  81.       RETURN
  82.       END
Advertisement
Add Comment
Please, Sign In to add comment