PROGRAM CO_SINTERING C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~C C Physical background : Co-sintering of two axi-symmetric particles C C.......................................................................C C 3-D Real Model written by Wendong Luo C C Date of start : 03/2012 C C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~C IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*11 JNAME CHARACTER*20 DFILE DIMENSION SPOINT(2),TIPCONST(4),XNS(1000,2),YNS(1000,2),DS(1000,2), + ANS(2,1000,2),AK(1000,2),FLUX(1000,2),VPORE(1000,2), + KEATEN(2),KADDN(2),BALL(2),AKMEAN(2) Pi=3.14159265359879323547d0 KBBENCH = 1 ! Set 'KBBENCH = 2' to run back-bench job CALL INPUT(JNAME,KJOB,KTBALL,ALMEANS0,SPOINT,Sigma_infy, + TIPCONST,XNS,YNS,ADDNCR,EATENCR, + NTSVO,NUMTS,NTSSURF,NTHIS,KBBENCH, + TIME,XC1,YC1,XC2,YC2) CALL OUTPUT(0,0,JNAME,KJOB,TIMESTEP,TIME, + ALMEANS0,SPOINT,TIPCONST,XNS,YNS) CALL MASS(KTBALL,SPOINT,XNS,YNS, ! in + BALL,AMASS) ! out PRINT *,' A_BALL(1) = ',BALL(1) PRINT *,' A_BALL(2) = ',BALL(2) BALL10=BALL(1) BALL20=BALL(2) DFILE=JNAME(:KJOB)//'.his' OPEN(21,FILE=DFILE) DFILE=JNAME(:KJOB)//'.eat' OPEN(24,FILE=DFILE) NOPT = NINT(FLOAT(NUMTS)/FLOAT(NTSVO)) NHIS = NINT(FLOAT(NUMTS)/FLOAT(NTHIS)) KOUTPUT = 0 IFILE = 0 FMASS_TIP= 0.0 FMASS_UP = 0.0 FMASS_LW = 0.0 AINT_AK_TIP = 0.0 KHIS = 0 WRITE(21,*) NTHIS-1 DO 1000 ITIME=1,NUMTS KOUTPUT=KOUTPUT+1 KHIS=KHIS+1 CALL ANALY(KTBALL,ALMEANS0,Sigma_infy,TIPCONST, ! in + SPOINT,XNS,YNS, ! in + DS,ANS,AK,FLUX,VPORE, ! out + UGRAIN,FLUX_TIP,T_FLUX_UP,T_FLUX_LW) ! out VPMAX=1.0D-10 DO 150 ISURF=1,2 NODES=IDNINT(SPOINT(ISURF)) AKMEAN(ISURF)=0.d0 DO 150 INODES=1,NODES IF(DABS(VPORE(INODES,ISURF)).GT.VPMAX) THEN VPMAX=DABS(VPORE(INODES,ISURF)) ENDIF IF(((ISURF.EQ.1).AND.(INODES.NE.NODES)).OR. + ((ISURF.EQ.2).AND.(INODES.NE.1))) THEN AKMEAN(ISURF)=AKMEAN(ISURF)+AK(INODES,ISURF)/(SPOINT(ISURF)-1.d0) ENDIF 150 CONTINUE TIMESTEP=ALMEANS0/(VPMAX*FLOAT(NTSSURF)) TIME=TIME+TIMESTEP FMASS_TIP= FMASS_TIP+FLUX_TIP*2.d0*Pi*XNS(1,2)*TIMESTEP FMASS_UP = FMASS_UP+T_FLUX_UP*2.d0*Pi*TIMESTEP FMASS_LW = FMASS_LW+T_FLUX_LW*2.d0*Pi*TIMESTEP AINT_AK_TIP = AINT_AK_TIP + AK(1,2)*TIMESTEP CALL UPDATE(UGRAIN,TIMESTEP,ALMEANS0,TIPCONST, ! in + VPORE,SPOINT,PORETIP,ANS, ! in + XNS,YNS,YC1) ! changed CALL REMESH(ITIME,ADDNCR,EATENCR,ALMEANS0,! in + TIPCONST,ANS, ! in + KEATEN,KADDN,SPOINT,XNS,YNS) ! changed CALL MASS(KTBALL,SPOINT,XNS,YNS, ! in + BALL,AMASS) ! out IF((ITIME.EQ.1).OR.(KHIS.EQ.NHIS)) THEN IF(KHIS.EQ.NHIS) KHIS=0 PRINT *,'ITIME = ',ITIME PRINT *,'TIME = ',TIME WRITE(*,160) UGRAIN,AK(1,2),XNS(1,2),DCENTER,AK(NODES-300,1), + AK(300,2),AKMEAN(1),AKMEAN(2) 160 FORMAT(2X,'UGRAIN= ',D25.15,2X,'K_tip= ',F20.10/ + 2X,'Neck Radius= ',F20.10,2X,'DCENTER= ',F20.10/ + 2X,'AK(NODES-300,1)= ',F25.15,2X,'AK(300,2)= ',F25.15/ + 2X,'K1_mean= ',F25.15,2X,'K2_mean= ',F25.15/) ERROR=DABS(AMASS-BALL10-BALL20)/(BALL10+BALL20) DCENTER=YC1-YC2 WRITE(21,*) ITIME WRITE(21,170) TIME WRITE(21,170) XNS(1,2),DCENTER WRITE(21,170) AK(1,2) !WRITE(21,170) AK(NODES-1,1),AK(NODES-2,1),AK(NODES-5,1) !WRITE(21,170) AK(2,2),AK(3,2),AK(6,2) !WRITE(21,170) AINT_AK_TIP,AKMEAN(1),AKMEAN(2) !WRITE(21,170) BALL(1),BALL(2) !WRITE(21,170) ERROR,AMASS !WRITE(21,180) FMASS_LW,FMASS_TIP,FMASS_UP !WRITE(21,180) VPORE(NODES-1,1),VPORE(NODES-2,1) !WRITE(21,180) VPORE(2,2),VPORE(3,2) WRITE(21,180) FLUX(NODES-1,1),FLUX(1,2) WRITE(21,180) FLUX(NODES-2,1),FLUX(2,2) WRITE(21,180) FLUX(NODES-3,1),FLUX(3,2) WRITE(21,180) FLUX(NODES-4,1),FLUX(4,2) WRITE(21,180) FLUX(NODES-5,1),FLUX(5,2) !WRITE(21,180) FLUX(NODES-10,1),FLUX(10,2) !WRITE(21,180) FLUX(NODES-15,1),FLUX(15,2) !WRITE(21,180) FLUX(NODES-20,1),FLUX(20,2) !WRITE(21,180) FLUX(NODES-25,1),FLUX(25,2) !WRITE(21,180) FLUX(NODES-30,1),FLUX(30,2) !WRITE(21,180) FLUX(NODES-50,1),FLUX(50,2) WRITE(21,180) UGRAIN 170 FORMAT(4F33.18) 180 FORMAT(4D33.20) ENDIF IF(KOUTPUT.EQ.NOPT) THEN IFILE=IFILE+1 KOUTPUT=0 CALL OUTPUT(ITIME,IFILE,JNAME,KJOB,TIMESTEP,TIME, + ALMEANS0,SPOINT,TIPCONST,XNS,YNS) ENDIF 1000 CONTINUE CLOSE(21) CLOSE(23) CLOSE(24) WRITE(*, 2000) 2000 FORMAT(/5X,'Normal termination.'/) STOP END SUBROUTINE MASS(KTBALL,SPOINT,XNS,YNS, ! in + BALL,AMASS) ! out IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),XNS(1000,2),YNS(1000,2),BALL(2) DO 500 IBALL=1,2 NODES=IDNINT(SPOINT(IBALL)) BALL(IBALL)=0.d0 DO 200 INS=1,NODES-1 R1 = XNS(INS,IBALL) R2 = XNS(INS+1,IBALL) AH = -YNS(INS+1,IBALL)+YNS(INS,IBALL) BALL(IBALL)=BALL(IBALL)+AH*(R1*R1+R1*R2+R2*R2) 200 CONTINUE BALL(IBALL)=BALL(IBALL)*3.14159265359879d0/3.d0 500 CONTINUE AMASS=BALL(1)+BALL(2) RETURN END SUBROUTINE INPUT(JNAME,KJOB,KTBALL,ALMEANS0,SPOINT,Sigma_infy, + TIPCONST,XNS,YNS,ADDNCR,EATENCR, + NTSVO,NUMTS,NTSSURF,NTHIS,KBBENCH, + TIME,XC1,YC1,XC2,YC2) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),TIPCONST(4),XNS(1000,2),YNS(1000,2) CHARACTER*11 JNAME,JNAMEUP CHARACTER*20 DFILE CHARACTER*4 CYES LOGICAL*4 EFILE 14 IF(KBBENCH.EQ.2) THEN OPEN(40,FILE='jobcard') READ(40,20) JNAME ELSE PRINT 15 15 FORMAT(/5X,'JOB NAME PLEASE :'/) READ(*,20) JNAME ENDIF 20 FORMAT(A11) KJOB=INDEX(JNAME,' ')-1 IF(KBBENCH.NE.2) THEN PRINT 45 45 FORMAT(5X,'Re-start ? ( y or n )'/) READ(*,50) CYES 50 FORMAT(A4) ELSE READ(40,50) CYES ENDIF IF(CYES(1:1).EQ.'Y'.OR.CYES(1:1).EQ.'y') THEN JNAME=JNAME(:KJOB)//'R' KJOB=KJOB+1 KRESTART=2 79 IF(KBBENCH.NE.2) THEN PRINT 80 80 FORMAT(/5X,'Restart data-file name please :'/) READ(*,90) JNAMEUP 90 FORMAT(A11) ELSE READ(40,90) JNAMEUP ENDIF KJOBR=INDEX(JNAMEUP,' ')-1 INQUIRE(FILE=JNAMEUP,EXIST=EFILE) IF(EFILE.NE..TRUE.) THEN PRINT 110, JNAMEUP 110 FORMAT(/5X,'Restart data file ',A11,' does not exist !'/ + /5X,'Did you type it incorrectly ? (y or n)'/) IF(KBBENCH.NE.2) THEN READ(*,50) CYES IF(CYES(1:1).EQ.'Y'.OR.CYES(1:1).EQ.'y') GOTO 79 PRINT 120 120 FORMAT(/5X,'Hi you fool ! what do you expect me to do then ?'/) ENDIF STOP ENDIF OPEN(17,FILE=JNAMEUP,STATUS='OLD') READ(17,*) ITIME,TIMESTEP,TIME READ(17,*) KTBALL DO 200 ISURF=1,2 READ(17,*) SPOINT(ISURF) DO 200 IPOINT=1,IDNINT(SPOINT(ISURF)) READ(17,*) JPOINT,XNS(IPOINT,ISURF),YNS(IPOINT,ISURF) 200 CONTINUE READ(17,*) ALMEANS0 CLOSE(17) ELSE KRESTART=0 TIME=0.0 DFILE=JNAME(:KJOB)//'.inp' INQUIRE(FILE=DFILE,EXIST=EFILE) IF(EFILE.NE..TRUE.) THEN PRINT 250, DFILE 250 FORMAT(/5X,'DATA FILE ', A20, + /5X,'DOES NOT EXIST ! WRONG JOB NAME ?' + /5X,'(y or n)'/) IF(KBBENCH.NE.2) THEN READ(*,50) CYES IF(CYES(1:1).EQ.'Y'.OR.CYES(1:1).EQ.'y') GOTO 14 PRINT 120 ENDIF STOP ENDIF OPEN(9,FILE=DFILE,STATUS='OLD') READ(9,*) XC1,YC1,XC2,YC2 READ(9,*) KTBALL,R1,R2,C DO 450 ISF=1,2 READ(9,*) SPOINT(ISF) DO 440 ISP=1,IDNINT(SPOINT(ISF)) READ(9,*) JPOINT,XNS(JPOINT,ISF),YNS(JPOINT,ISF) 440 CONTINUE 450 CONTINUE READ(9,*) ALMEANS0 CLOSE(9) ENDIF DFILE=JNAME(:KJOB)//'.run' OPEN(23,FILE=DFILE) WRITE(23,460) 460 FORMAT(/12X,'******** RUNNING TIME INFORMATION ********'//) NUMTS=10000000 NTSSURF=3 TIPCONST(1)=6.8531648D-34 !!!DDs TIPCONST(4)=6.8531648D-36 !!!DDb TIPCONST(2)=1.1d0 !!!Gamma_s Fai=3.1415926535898d0*0.3611111111111111 !Cos(FAI)=Gamma_gb/2Gamma_s,65 TIPCONST(3)=2.D0*DCOS(Fai) EATENCR = 0.05 ADDNCR = 1.2 NTSVO = 100 Sigma_infy = 0.0 IF(NUMTS.LE.10000) THEN NTHIS = NUMTS ELSE NTHIS = 1000 ENDIF WRITE(*, 700) JNAME,KTBALL WRITE(23,700) JNAME,KTBALL 700 FORMAT(/5X,'Job name : ',A11,' KTBALL = ',I3/) IF(KRESTART.EQ.2) THEN WRITE(*, 750) JNAMEUP WRITE(23,750) JNAMEUP 750 FORMAT(/5X,'This is a restart job from data file : ',A11/) ENDIF 760 WRITE(* ,770) NUMTS,NTSSURF,(TIPCONST(JTIP),JTIP=1,4), + Sigma_infy,EATENCR,ADDNCR,NTHIS WRITE(23,770) NUMTS,NTSSURF,(TIPCONST(JTIP),JTIP=1,4), + Sigma_infy,EATENCR,ADDNCR,NTHIS 770 FORMAT(/5X,'1. Total time step numbers : ',I10 + /5X,'2. Time-step length control : ',I10 + /5X,'3. DDs : ',D20.10 + /5X,'4. Gamma_s : ',F20.10 + /5X,'5. Gamma_b/Gamma_s : ',F20.10 + /5X,'6. DDb : ',D20.10 + /5X,'7. Remote Stress : ',D20.10 + /5X,'8. Minimum surface-element length : ',F20.10 + /5X,'9. Maximum surface-element length : ',F20.10 + /5X,'10. Number of writtings to *.his file : ',I6/) IF(KBBENCH.NE.2) THEN WRITE(*,780) 780 FORMAT(/5X,'Do you want to change them ? (y or n )'/) READ(*,50) CYES IF(CYES(1:1).EQ.'Y'.OR.CYES(1:1).EQ.'y') THEN WRITE(*,790) 790 FORMAT(/5X,'Which one ?'/) READ(*,*) ICHANGE WRITE(*,800) 800 FORMAT(/5X,'Input the new value :'/) IF(ICHANGE.EQ.1) READ(*,*) NUMTS IF(ICHANGE.EQ.2) READ(*,*) NTSSURF IF(ICHANGE.EQ.3) READ(*,*) TIPCONST(1) IF(ICHANGE.EQ.4) READ(*,*) TIPCONST(2) IF(ICHANGE.EQ.5) READ(*,*) TIPCONST(3) IF(ICHANGE.EQ.6) READ(*,*) TIPCONST(4) IF(ICHANGE.EQ.7) READ(*,*) Sigma_infy IF(ICHANGE.EQ.8) READ(*,*) EATENCR IF(ICHANGE.EQ.9) READ(*,*) ADDNCR IF(ICHANGE.EQ.10) READ(*,*) NTHIS GOTO 760 ENDIF ENDIF WRITE(23,870) R1,R2,C 870 FORMAT(/5X,'11. Radius of upper ball : ',F30.20/ + 5X,'12. Radius of lower ball : ',F30.20/ + 5X,'13. Initial Neck size : ',F30.20) RETURN END SUBROUTINE OUTPUT(ITIME,IFILE,JNAME,KJOB,TIMESTEP,TIME, + ALMEANS0,SPOINT,TIPCONST,XNS,YNS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION TIPCONST(4),SPOINT(2),XNS(1000,2),YNS(1000,2) CHARACTER*11 JNAME CHARACTER*20 DFILE CALL FILENAME(IFILE,JNAME,KJOB,DFILE) OPEN(17,FILE=DFILE) WRITE(17,*) ITIME,TIMESTEP,TIME C WRITE(17,*) KTBALL DO 200 ISURF=1,2 WRITE(17,*) SPOINT(ISURF) DO 200 IPOINT=1,IDNINT(SPOINT(ISURF)) WRITE(17,*) IPOINT,XNS(IPOINT,ISURF),YNS(IPOINT,ISURF) 200 CONTINUE C WRITE(17,*) ALMEANS0 CLOSE(17) C300 FORMAT(F20.10,TR6,F20.10) RETURN END C*******************************************************************C C UPDATING BLOCK C C*******************************************************************C SUBROUTINE UPDATE(UGRAIN,TIMESTEP,ALMEANS0,TIPCONST, ! in + VPORE,SPOINT,PORETIP,ANS, ! in + XNS,YNS,YC1) ! changed IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),TIPCONST(4),VPORE(1000,2), + XNS(1000,2),YNS(1000,2),ANS(2,1000,2) CALL UPPORE(TIMESTEP,SPOINT,ANS,VPORE, ! in + XNS,YNS) ! changed CALL NEWXY(TIMESTEP,UGRAIN,ALMEANS0,SPOINT,TIPCONST,ANS, ! in + XNS,YNS,YC1) ! changed C YSYM=YNS(1,2) C NODES=IDNINT(SPOINT(1)) C DO 500 INS=1,NODES-1 C XNS(NODES-INS,1)=XNS(INS+1,2) C YNS(NODES-INS,1)=YSYM+(YSYM-YNS(INS+1,2)) C500 CONTINUE RETURN END SUBROUTINE UPPORE(TIMESTEP,SPOINT,ANS,VPORE, ! in + XNS,YNS) ! changed IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),ANS(2,1000,2),VPORE(1000,2), + XNS(1000,2),YNS(1000,2) DO 900 ISURF=1,2 NODES=IDNINT(SPOINT(ISURF)) IF(ISURF.EQ.1) THEN DO 800 INODES=1,NODES-1 XNS(INODES,ISURF)=XNS(INODES,ISURF)+ + TIMESTEP*ANS(1,INODES,ISURF) + *VPORE(INODES,ISURF) YNS(INODES,ISURF)=YNS(INODES,ISURF)+ + TIMESTEP*ANS(2,INODES,ISURF) + *VPORE(INODES,ISURF) 800 CONTINUE ELSE DO 850 INODES=2,NODES XNS(INODES,ISURF)=XNS(INODES,ISURF)+ + TIMESTEP*ANS(1,INODES,ISURF) + *VPORE(INODES,ISURF) YNS(INODES,ISURF)=YNS(INODES,ISURF)+ + TIMESTEP*ANS(2,INODES,ISURF) + *VPORE(INODES,ISURF) 850 CONTINUE ENDIF 900 CONTINUE RETURN END SUBROUTINE NEWXY(TIMESTEP,UGRAIN,ALMEANS0,SPOINT,TIPCONST,ANS, ! in + XNS,YNS,YC1) ! changed IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),TIPCONST(4), + XNS(1000,2),YNS(1000,2),ANS(2,1000,2), + XLIN1(2),YLIN1(2),XLIN2(2),YLIN2(2), + GBX(2),GBY(2) NODES1=IDNINT(SPOINT(1)) DO 100 INS=1,NODES1-1 100 YNS(INS,1)=YNS(INS,1)+TIMESTEP*UGRAIN YC1=YC1+TIMESTEP*UGRAIN C- Find the new tip position DIST=100.d0*ALMEANS0 FAI=DACOS(TIPCONST(3)/2.d0) GBX(1)=0.d0 GBX(2)=XNS(1,2) GBY(1)=YNS(1,2) GBY(2)=GBY(1) CALL DRAWLINE(DIST,FAI,GBX,GBY, + XNS(NODES1-1,1),YNS(NODES1-1,1), + XNS(2,2),YNS(2,2), + XLIN1,YLIN1,XLIN2,YLIN2) CALL INTERSECT(XLIN1,YLIN1,XLIN2,YLIN2, + XJOINT,YJOINT,NOJOINT) IF(NOJOINT.NE.2) THEN XNS(NODES1,1)=XJOINT YNS(NODES1,1)=YJOINT XNS(1,2)=XJOINT YNS(1,2)=YJOINT ELSE PRINT *,'Error A in UPDATING' STOP ENDIF RETURN END SUBROUTINE REMESH(ITIME,ADDNCR,EATENCR,ALMEANS0,! in + TIPCONST,ANS, ! in + KEATEN,KADDN,SPOINT,XNS,YNS) ! changed IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),XNS(1000,2),YNS(1000,2),ANS(2,1000,2), + TIPCONST(4),KEATEN(2),XLIN1(2),YLIN1(2),XLIN2(2), + YLIN2(2),GBX(2),GBY(2),XSURF3(2),YSURF3(2), + KADDN(2) FAI=DACOS(TIPCONST(3)/2.d0) COSFAI=TIPCONST(3)/2.d0 SINFAI=DSQRT(1.d0-COSFAI*COSFAI) C----------- add node if a surface element becomes too long DO 500 ISURF=1,2 KADDN(ISURF)=0 50 NODES=IDNINT(SPOINT(ISURF)) DO 490 IES=1,NODES-1 ALIES=DSQRT((XNS(IES,ISURF)-XNS(IES+1,ISURF))**2 + +(YNS(IES,ISURF)-YNS(IES+1,ISURF))**2) IF((ALIES/ALMEANS0).GT.ADDNCR) THEN KADDN(ISURF)=KADDN(ISURF)+1 SPOINT(ISURF)=SPOINT(ISURF)+1.d0 IF(SPOINT(ISURF).GT.999.d0) THEN PRINT *,'Surface nodes > 1000 !' STOP ENDIF XNEW=(XNS(IES,ISURF)+XNS(IES+1,ISURF))/2.d0 YNEW=(YNS(IES,ISURF)+YNS(IES+1,ISURF))/2.d0 DO 400 INS=NODES,IES+1,-1 XNS(INS+1,ISURF)=XNS(INS,ISURF) 400 YNS(INS+1,ISURF)=YNS(INS,ISURF) XNS(IES+1,ISURF)=XNEW YNS(IES+1,ISURF)=YNEW GOTO 50 ENDIF 490 CONTINUE 500 CONTINUE C------- delete element which is shorter than prescribed length ----C DO 3000 ISURF=1,2 KEATEN(ISURF)=0 NODES=IDNINT(SPOINT(ISURF)) DO 2000 IES=1,NODES-1 XSU1=XNS(IES,ISURF) XSU2=XNS(IES+1,ISURF) YSU1=YNS(IES,ISURF) YSU2=YNS(IES+1,ISURF) ALSU=DSQRT((XSU1-XSU2)**2+(YSU1-YSU2)**2) IF((ALSU/ALMEANS0).LT.EATENCR) THEN KEATEN(ISURF)=KEATEN(ISURF)+1 IF((ISURF.EQ.2).AND.(IES.EQ.1)) THEN GBX(1)=0.d0 GBX(2)=XNS(1,2) GBY(1)=YNS(1,2) GBY(2)=YNS(1,2) CALL DRAWLINE(100.d0*ALMEANS0,FAI,GBX,GBY, + XNS(1,2),YNS(1,2),XNS(1,2),YNS(1,2), + XLIN1,YLIN1,XLIN2,YLIN2) XSURF3(1)=XNS(3,ISURF) XSURF3(2)=XNS(3,ISURF)+ANS(1,3,ISURF)*100.d0*ALMEANS0 YSURF3(1)=YNS(3,ISURF) YSURF3(2)=YNS(3,ISURF)+ANS(2,3,ISURF)*100.d0*ALMEANS0 CALL INTERSECT(XLIN2,YLIN2, + XSURF3,YSURF3, + XJOINTN,YJOINTN,NOJOINTN) IF(NOJOINTN.NE.2) THEN TEMP1=TAREA(XNS(2,ISURF),YNS(2,ISURF), + XNS(3,ISURF),YNS(3,ISURF), + XJOINTN,YJOINTN) TEMP2=TAREA(XNS(4,ISURF),YNS(4,ISURF), + XNS(3,ISURF),YNS(3,ISURF), + XJOINTN,YJOINTN) ERRORN = TEMP1+TEMP2 ELSE ERRORN = 1.0D+20 ENDIF XSURF3(1)=XNS(3,ISURF) XSURF3(2)=XNS(4,ISURF) YSURF3(1)=YNS(3,ISURF) YSURF3(2)=YNS(4,ISURF) CALL INTERSECT(XLIN2,YLIN2, + XSURF3,YSURF3, + XJOINTT,YJOINTT,NOJOINTT) IF(NOJOINTT.NE.2) THEN ERRORT=TAREA(XNS(2,ISURF),YNS(2,ISURF), + XNS(3,ISURF),YNS(3,ISURF), + XJOINTT,YJOINTT) ELSE ERRORT = 1.0D+20 ENDIF IF((NOJOINTN.NE.2).OR.(NOJOINTT.NE.2)) THEN IF(ERRORT.GT.ERRORN) THEN XNS(3,ISURF)=XJOINTN YNS(3,ISURF)=YJOINTN ELSE XNS(3,ISURF)=XJOINTT YNS(3,ISURF)=YJOINTT ENDIF ELSE PRINT *,'ERROR Q1 in eatnodes' STOP ENDIF XNS(2,ISURF)=DABS(ALMEANS0)*1.0D+8 YNS(2,ISURF)=DABS(ALMEANS0)*1.0D+8 ELSE IF((ISURF.EQ.1).AND.(IES.EQ.NODES-1)) THEN GBX(1)=0.d0 GBX(2)=XNS(1,2) GBY(1)=YNS(1,2) GBY(2)=YNS(1,2) CALL DRAWLINE(100.d0*ALMEANS0,FAI,GBX,GBY, + XNS(1,2),YNS(1,2),XNS(1,2),YNS(1,2), + XLIN1,YLIN1,XLIN2,YLIN2) XSURF3(1)=XNS(NODES-2,ISURF) XSURF3(2)=XNS(NODES-2,ISURF)+ + ANS(1,NODES-2,ISURF)*100.d0*ALMEANS0 YSURF3(1)=YNS(NODES-2,ISURF) YSURF3(2)=YNS(NODES-2,ISURF)+ + ANS(2,NODES-2,ISURF)*100.d0*ALMEANS0 CALL INTERSECT(XLIN1,YLIN1, + XSURF3,YSURF3, + XJOINTN,YJOINTN,NOJOINTN) IF(NOJOINTN.NE.2) THEN TEMP1=TAREA(XNS(NODES-1,ISURF),YNS(NODES-1,ISURF), + XNS(NODES-2,ISURF),YNS(NODES-2,ISURF), + XJOINTN,YJOINTN) TEMP2=TAREA(XNS(NODES-2,ISURF),YNS(NODES-2,ISURF), + XNS(NODES-3,ISURF),YNS(NODES-3,ISURF), + XJOINTN,YJOINTN) ERRORN = TEMP1+TEMP2 ELSE ERRORN = 1.0D+20 ENDIF XSURF3(1)=XNS(NODES-2,ISURF) XSURF3(2)=XNS(NODES-3,ISURF) YSURF3(1)=YNS(NODES-2,ISURF) YSURF3(2)=YNS(NODES-3,ISURF) CALL INTERSECT(XLIN1,YLIN1, + XSURF3,YSURF3, + XJOINTT,YJOINTT,NOJOINTT) IF(NOJOINTT.NE.2) THEN ERRORT=TAREA(XNS(NODES-1,ISURF),YNS(NODES-1,ISURF), + XNS(NODES-2,ISURF),YNS(NODES-2,ISURF), + XJOINTT,YJOINTT) ELSE ERRORT = 1.0D+20 ENDIF IF((NOJOINTN.NE.2).OR.(NOJOINTT.NE.2)) THEN IF(ERRORT.GT.ERRORN) THEN XNS(NODES-2,ISURF)=XJOINTN YNS(NODES-2,ISURF)=YJOINTN ELSE XNS(NODES-2,ISURF)=XJOINTT YNS(NODES-2,ISURF)=YJOINTT ENDIF ELSE PRINT *,'ERROR Q2 in eatnodes' STOP ENDIF XNS(NODES-1,ISURF)=DABS(ALMEANS0)*1.0D+8 YNS(NODES-1,ISURF)=DABS(ALMEANS0)*1.0D+8 ELSE XNS(IES+1,ISURF)=(XNS(IES,ISURF)+XNS(IES+1,ISURF))/2.d0 YNS(IES+1,ISURF)=(YNS(IES,ISURF)+YNS(IES+1,ISURF))/2.d0 XNS(IES,ISURF)=DABS(ALMEANS0)*1.0D+8 YNS(IES,ISURF)=DABS(ALMEANS0)*1.0D+8 ENDIF ENDIF 2000 CONTINUE IF(KEATEN(ISURF).GT.0) THEN DO 2500 INODES=1,NODES IF(XNS(INODES,ISURF).GT.DABS(ALMEANS0)*1.0D+7) THEN IF(INODES.LT.NODES) THEN DO 2400 JNODES=INODES,NODES-1 XNS(JNODES,ISURF)=XNS(JNODES+1,ISURF) 2400 YNS(JNODES,ISURF)=YNS(JNODES+1,ISURF) ENDIF ENDIF 2500 CONTINUE SPOINT(ISURF)=SPOINT(ISURF)-FLOAT(KEATEN(ISURF)) ENDIF 3000 CONTINUE IF(KEATEN(1).GT.0) THEN PRINT *,' ** KEATEN(1) = ',KEATEN(1),' **' WRITE(24,651) ITIME,KEATEN(1) 651 FORMAT(5X,'ITIME = ',I10,' KEATEN(1) = ',I5) ENDIF IF(KEATEN(2).GT.0) THEN PRINT *,' ** KEATEN(2) = ',KEATEN(2),' **' WRITE(24,652) ITIME,KEATEN(2) 652 FORMAT(5X,'ITIME = ',I10,' KEATEN(2) = ',I5) ENDIF IF(KADDN(1).GT.0) THEN PRINT *,' ** KADDN(1) = ',KADDN(1),' **' WRITE(24,653) ITIME,KADDN(1) 653 FORMAT(5X,'ITIME = ',I10,' KADDN(1) = ',I5) ENDIF IF(KADDN(2).GT.0) THEN PRINT *,' ** KADDN(2) = ',KADDN(2),' **' WRITE(24,654) ITIME,KADDN(2) 654 FORMAT(5X,'ITIME = ',I10,' KADDN(2) = ',I5) ENDIF RETURN END FUNCTION TAREA(X1,Y1,X2,Y2,X3,Y3) IMPLICIT DOUBLE PRECISION (A-H,O-Z) AL12 = DSQRT((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1)) AL13 = DSQRT((X3-X1)*(X3-X1)+(Y3-Y1)*(Y3-Y1)) AL23 = DSQRT((X2-X3)*(X2-X3)+(Y2-Y3)*(Y2-Y3)) P1 = 0.5D0*(AL12+AL23+AL13) TAREA = DSQRT(P1*(P1-AL12)*(P1-AL23)*(P1-AL13)) RETURN END SUBROUTINE DRAWLINE(DIST,FAI,XLIN0,YLIN0, + X1,Y1,X2,Y2, + XLIN1,YLIN1,XLIN2,YLIN2) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XLIN0(2),YLIN0(2),XLIN1(2),YLIN1(2), + XLIN2(2),YLIN2(2) XLIN1(1)=X1 YLIN1(1)=Y1 XLIN2(1)=X2 YLIN2(1)=Y2 AL0=DSQRT((XLIN0(1)-XLIN0(2))**2+(YLIN0(1)-YLIN0(2))**2) ALPHA=DACOS(DABS(XLIN0(1)-XLIN0(2))/AL0) WHICHLIN0=(XLIN0(1)-XLIN0(2))*(YLIN0(1)-YLIN0(2)) IF(WHICHLIN0.GE.0) THEN BETA1=FAI+ALPHA BETA2=FAI-ALPHA ELSE BETA1=FAI-ALPHA BETA2=FAI+ALPHA ENDIF XLIN1(2)=X1+DIST*DCOS(BETA1) YLIN1(2)=Y1+DIST*DSIN(BETA1) XLIN2(2)=X2+DIST*DCOS(BETA2) YLIN2(2)=Y2-DIST*DSIN(BETA2) RETURN END SUBROUTINE INTERSECT(XLINE1,YLINE1,XLINE2,YLINE2, + XJOINT,YJOINT,NOJOINT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XLINE1(2),YLINE1(2),XLINE2(2),YLINE2(2) AL1=DSQRT((XLINE1(1)-XLINE1(2))**2+ + (YLINE1(1)-YLINE1(2))**2) AL2=DSQRT((XLINE2(1)-XLINE2(2))**2+ + (YLINE2(1)-YLINE2(2))**2) VERT1=DABS(XLINE1(1)-XLINE1(2))/AL1 VERT2=DABS(XLINE2(1)-XLINE2(2))/AL2 NOJOINT=0 IF((VERT1.LT.1.0D-4).AND.(VERT2.LT.1.0D-4)) THEN NOJOINT=2 RETURN ELSE IF(VERT1.LT.1.0D-4) THEN L1VERT=2 L2VERT=0 AKL2=(YLINE2(2)-YLINE2(1))/(XLINE2(2)-XLINE2(1)) ELSE IF(VERT2.LT.1.0D-4) THEN L1VERT=0 L2VERT=2 AKL1=(YLINE1(2)-YLINE1(1))/(XLINE1(2)-XLINE1(1)) ELSE L1VERT=0 L2VERT=0 AKL1=(YLINE1(2)-YLINE1(1))/(XLINE1(2)-XLINE1(1)) AKL2=(YLINE2(2)-YLINE2(1))/(XLINE2(2)-XLINE2(1)) ENDIF IF((L1VERT.EQ.0).AND.(L2VERT.EQ.0)) THEN IF(DABS(AKL1-AKL2).LT.1.0D-3) THEN NOJOINT=2 RETURN ENDIF ENDIF IF(L1VERT.EQ.2) THEN XJOINT=XLINE1(1) YJOINT=AKL2*(XJOINT-XLINE2(1))+YLINE2(1) ELSE IF(L2VERT.EQ.2) THEN XJOINT=XLINE2(1) YJOINT=AKL1*(XJOINT-XLINE1(1))+YLINE1(1) ELSE BL1=YLINE1(1)-AKL1*XLINE1(1) BL2=YLINE2(1)-AKL2*XLINE2(1) XJOINT=-(BL1-BL2)/(AKL1-AKL2) YJOINT=AKL1*XJOINT+BL1 ENDIF RETURN END C**************************************************************************C C Analysis Block C C**************************************************************************C SUBROUTINE ANALY(KTBALL,ALMEANS0,Sigma_infy,TIPCONST, ! in + SPOINT,XNS,YNS, ! in + DS,ANS,AK,FLUX,VPORE, ! out + UGRAIN,FLUX_TIP,T_FLUX_UP,T_FLUX_LW) ! out IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),TIPCONST(4),XNS(1000,2),YNS(1000,2),DS(1000,2), + ANS(2,1000,2),AK(1000,2),FLUX(1000,2),VPORE(1000,2) C CHARACTER*4 CYES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C TIPCONST(1) : DDb C C TIPCONST(2) : Gamma_s C C TIPCONST(3) : Gamma_b/Gamma_s C C TIPCONST(4) : DDb C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL CLEAR(DS,2*1000) CALL CLEAR(ANS,4*1000) CALL CLEAR(AK,2*1000) CALL CLEAR(FLUX,2*1000) CALL CLEAR(VPORE,2*1000) ! Cs = TIPCONST(1)/TIPCONST(4)/100 DDs = TIPCONST(1) DDb = TIPCONST(4) Gamma_s = TIPCONST(2) COSFAi = TIPCONST(3)/2.d0 SINFAI = DSQRT(1.d0-COSFAI**2) CALL COMP_DS(SPOINT,XNS,YNS, ! in + DS) ! out CALL COMP_AK_ANS(KTBALL,ALMEANS0,SPOINT,XNS,YNS, ! in + ANS,AK) ! out CALL COMP_FLUX(DDs,SPOINT,AK,DS, ! in + FLUX) ! out CALL COMP_VPORE(SPOINT,XNS,DS,FLUX, ! in + VPORE) ! out CALL COMPLETE(DDs,DDb,Gamma_s,SINFAI,Sigma_infy, ! in + SPOINT,XNS,DS, ! in + AK,FLUX,VPORE) ! complete R_NECK=XNS(1,2) PI=3.1415926535897932d0 Sigma_C=AK(1,2)*Gamma_s UGRAIN = DDb*8.d0*((Sigma_infy/R_NECK**2)-(Sigma_C/R_NECK**2)- + (2.d0*Gamma_s*SINFAI/R_NECK**3)) FLUX_TIP = -0.5d0*UGRAIN*R_NECK NODES=IDNINT(SPOINT(1)) T_FLUX_UP=FLUX(NODES-1,1)*(XNS(NODES-1,1)+XNS(NODES,1))/2.d0 T_FLUX_LW=FLUX(1,2)*(XNS(1,2)+XNS(2,2))/2.d0 C TLEFT=TEMP1+TEMP2+TEMP3 C TEMP4=VPORE(2,2)*((XNS(1,2)+0.5d0*(XNS(1,2)+XNS(2,2)))/2.d0) C + *DS(1,2)/2.d0 C TEMP5=VPORE(NODES-1,1)*((XNS(NODES,1)+0.5d0*(XNS(NODES-1,1)+ C + XNS(NODES,1)))/2.d0)*DS(NODES-1,1)/2.d0 C TRIGHT=TEMP4 + TEMP5 C PRINT *,'TLEFT = ',TLEFT C PRINT *,'TRIGHT = ',TRIGHT C PRINT *,'FLUX_TIP, FLUX(1,2), FLUX(NODES-1,1) = ' c PRINT *, FLUX_TIP, FLUX(1,2), FLUX(NODES-1,1) C PRINT *,'VPORE(2,2), VPORE(NODES-1,1)' C PRINT *, VPORE(2,2), VPORE(NODES-1,1) C READ(*,50) CYES C50 FORMAT(A4) RETURN END SUBROUTINE COMP_DS(SPOINT,XNS,YNS, ! in + DS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),XNS(1000,2),YNS(1000,2),DS(1000,2) DO 1000 ISURF=1,2 NODES=IDNINT(SPOINT(ISURF)) DO 900 INS=1,NODES-1 DS(INS,ISURF)=DSQRT((XNS(INS,ISURF)-XNS(INS+1,ISURF))**2+ + (YNS(INS,ISURF)-YNS(INS+1,ISURF))**2) 900 CONTINUE 1000 CONTINUE RETURN END SUBROUTINE COMP_AK_ANS(KTBALL,ALMEANS0,SPOINT,XNS,YNS, ! in + ANS,AK) ! out IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),XNS(1000,2),YNS(1000,2), + ANS(2,1000,2),AK(1000,2) DO 1000 ISURF=1,2 NODES=IDNINT(SPOINT(ISURF)) IF(ISURF.EQ.1) THEN NODESF=1 NODESL=NODES-1 ELSE IF(ISURF.EQ.2) THEN NODESF=2 NODESL=NODES ENDIF DO 1000 IN=NODESF,NODESL IF((ISURF.EQ.1).AND.(IN.EQ.1)) THEN IF(KTBALL.EQ.2) THEN X1 = XNS(IN,ISURF) Y1 = YNS(IN,ISURF) XL = -XNS(IN+1,ISURF) YL = +YNS(IN+1,ISURF) X2 = XNS(IN+1,ISURF) Y2 = YNS(IN+1,ISURF) ELSE X1 = XNS(IN,ISURF) Y1 = YNS(IN,ISURF) XL = XNS(IN+1,ISURF) YL = YNS(IN,ISURF)+(YNS(IN,ISURF)-YNS(IN+1,ISURF)) X2 = XNS(IN+1,ISURF) Y2 = YNS(IN+1,ISURF) ENDIF ELSE IF((ISURF.EQ.2).AND.(IN.EQ.NODES)) THEN IF(KTBALL.EQ.2) THEN X1 = XNS(IN,ISURF) Y1 = YNS(IN,ISURF) XL = XNS(IN-1,ISURF) YL = YNS(IN-1,ISURF) X2 = -XNS(IN-1,ISURF) Y2 = +YNS(IN-1,ISURF) ELSE X1 = XNS(IN,ISURF) Y1 = YNS(IN,ISURF) XL = XNS(IN-1,ISURF) YL = YNS(IN-1,ISURF) X2 = XNS(IN-1,ISURF) Y2 = YNS(IN,ISURF)-(YNS(IN-1,ISURF)-YNS(IN,ISURF)) ENDIF ELSE IF((IN.NE.1).AND.(IN.NE.NODES)) THEN X1 = XNS(IN,ISURF) Y1 = YNS(IN,ISURF) XL = XNS(IN-1,ISURF) YL = YNS(IN-1,ISURF) X2 = XNS(IN+1,ISURF) Y2 = YNS(IN+1,ISURF) ENDIF CALL CAK1(X1,Y1,XL,YL,X2,Y2,ALMEANS0, ! in + ANS1,ANS2,AK1) ! out ANS(1,IN,ISURF) = ANS1 ANS(2,IN,ISURF) = ANS2 IF((DABS(XNS(IN,ISURF))/ALMEANS0).GT.1.0D-8) THEN AK2 = -ANS1/XNS(IN,ISURF) ELSE AK2 = AK1 ENDIF AK(IN,ISURF) = AK1 + AK2 1000 CONTINUE RETURN END SUBROUTINE CAK1(X1,Y1,XL,YL,X2,Y2,ALMEANS0, ! in + ANS1,ANS2,AK1) ! out IMPLICIT DOUBLE PRECISION (A-H,O-Z) AL = DSQRT((Y2-YL)**2+(X2-XL)**2) ATx = (X2-XL)/AL ATy = (Y2-YL)/AL ANx = -ATy ANy = ATx ALINE = (YL-Y1)*(X2-X1)-(Y2-Y1)*(XL-X1) Stra = DSQRT(DABS(ALINE))/ALMEANS0 IF(Stra.LT.1.0D-7) THEN ! Nodes L,1,2 form a straight line ANS1 = ANx ANS2 = ANy QL = 0.d0 ELSE TEMP1 = XL*XL-X1*X1+YL*YL-Y1*Y1 TEMP2 = X2*X2-X1*X1+Y2*Y2-Y1*Y1 Xr = -0.5D0*((Y2-Y1)*TEMP1-(YL-Y1)*TEMP2)/ALINE Yr = 0.5D0*((X2-X1)*TEMP1-(XL-X1)*TEMP2)/ALINE R = DSQRT((Yr-Y1)**2+(Xr-X1)**2) AN1x0 = (Xr-X1)/R AN1y0 = (Yr-Y1)/R CONCAVE = AN1x0*ANx + AN1y0*ANy IF(CONCAVE.GT.0.D0) THEN ! Concave surface at N1 ANS1 = AN1x0 ANS2 = AN1y0 AK1 = 1.d0/R ELSE ! Convax surface at N1 ANS1 = -AN1x0 ANS2 = -AN1y0 AK1 = -1.d0/R ENDIF ENDIF 4900 CONTINUE RETURN END SUBROUTINE COMP_FLUX(DDs,SPOINT,AK,DS, ! in + FLUX) ! out IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),DS(1000,2),AK(1000,2),FLUX(1000,2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Has problem DO 400 ISURF=1,2 NODES=IDNINT(SPOINT(ISURF)) IF(ISURF.EQ.1) THEN NODESF=1 NODESL=NODES-2 ELSE IF(ISURF.EQ.2) THEN NODESF=2 NODESL=NODES-1 ENDIF DO 200 IN = NODESF,NODESL FLUX(IN,ISURF) = DDs*(AK(IN+1,ISURF)-AK(IN,ISURF))/DS(IN,ISURF) 200 CONTINUE 400 CONTINUE RETURN END !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE COMP_VPORE(SPOINT,XNS,DS,FLUX, ! in + VPORE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),XNS(1000,2),DS(1000,2),FLUX(1000,2), + VPORE(1000,2) DO 1000 ISURF=1,2 NODES=IDNINT(SPOINT(ISURF)) IF(ISURF.EQ.1) THEN NODESF=2 NODESL=NODES-2 ELSE IF(ISURF.EQ.2) THEN NODESF=3 NODESL=NODES-1 ENDIF DO 800 IN=NODESF,NODESL TEMP1 = XNS(IN-1,ISURF)+XNS(IN,ISURF) TEMP2 = XNS(IN+1,ISURF)+XNS(IN,ISURF) TEMP3 = XNS(IN-1,ISURF)+XNS(IN,ISURF)*2.D0+XNS(IN+1,ISURF) TEMP4 = DS(IN-1,ISURF)+DS(IN,ISURF) VPORE(IN,ISURF) = 4.D0*(TEMP1*FLUX(IN-1,ISURF)- + TEMP2*FLUX(IN,ISURF)) + /(TEMP3*TEMP4) 800 CONTINUE IF(ISURF.EQ.1) THEN VPORE(1,ISURF) = -4.d0*FLUX(1,ISURF)/DS(1,ISURF) ELSE IF(ISURF.EQ.2) THEN NODES=IDNINT(SPOINT(ISURF)) VPORE(NODES,ISURF) = 4.d0*FLUX(NODES-1,ISURF)/ + DS(NODES-1,ISURF) ENDIF 1000 CONTINUE RETURN END SUBROUTINE COMPLETE(DDs,DDb,Gamma_s,SINFAI,Sigma_infy, ! in !lack of Cs + SPOINT,XNS,DS, ! in + AK,FLUX,VPORE) ! complete IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SPOINT(2),XNS(1000,2),DS(1000,2),FLUX(1000,2), + VPORE(1000,2),AK(1000,2) NODES = IDNINT(SPOINT(1)) ! F1_UP = 4.d0*((XNS(1,2)+XNS(2,2))*Cs*AK(2,2)/DS(1,2) ! + -(XNS(2,2)+XNS(3,2))*FLUX(2,2)) ! ! F1_BT = (XNS(1,2)+2.d0*XNS(2,2)+XNS(3,2))*(DS(1,2)+DS(2,2)) ! ! F1 = F1_UP/F1_BT ! ! F2 = -4.d0*(XNS(1,2)+XNS(2,2))*Cs/(DS(1,2)*F1_BT) ! ! F3_UP = 4.d0*((XNS(NODES-1,1)+XNS(NODES,1))*Cs*AK(NODES-1,1)/ ! + DS(NODES-1,1) ! + +(XNS(NODES-2,1)+XNS(NODES-1,1))*FLUX(NODES-2,1)) ! ! F3_BT = (XNS(NODES-2,1)+2.d0*XNS(NODES-1,1)+XNS(NODES,1))* ! + (DS(NODES-2,1)+DS(NODES-1,1)) ! ! F3 = F3_UP/F3_BT ! ! F4 = -4.d0*(XNS(NODES-1,1)+XNS(NODES,1))*Cs/ ! + (DS(NODES-1,1)*F3_BT) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! A1 = DDs/DS(NODES-1,1)!DDs*(XNS(NODES-1,1)+XNS(NODES,1))/(2.d0*DS(NODES-1,1)) A2 = DDs/DS(1,2)!DDs*(XNS(1,2)+XNS(2,2))/(2.d0*DS(1,2)) A3 = 4.d0*DDb*Gamma_s/XNS(1,2) ! A4 = -F2*(3.d0*XNS(1,2)+XNS(2,2))*DS(1,2)/8.d0 ! A5 = -F4*(3.d0*XNS(NODES,1)+XNS(NODES-1,1))*DS(NODES-1,1)/8.d0 A = A1 + A2 + A3 ! B1 = F1*(3.d0*XNS(1,2)+XNS(2,2))*DS(1,2)/8.d0 ! B2 = F3*(3.d0*XNS(NODES,1)+XNS(NODES-1,1))*DS(NODES-1,1)/8.d0 B3 = DDs*AK(NODES-1,1)/DS(NODES-1,1)!A1*AK(NODES-1,1) B4 = DDs*AK(2,2)/DS(1,2) !A2*AK(2,2) B5 = 4.d0*DDb*Sigma_infy/XNS(1,2) B6 = -8.d0*DDb*Gamma_s*SINFAI/(XNS(1,2)**2) B = B3 + B4 + B5 + B6 AK(NODES,1) = B/A AK(1,2) = B/A FLUX(NODES-1,1) = DDs*(AK(NODES,1)-AK(NODES-1,1))/DS(NODES-1,1)!!!same IN = NODES-1 ISURF = 1 TEMP1 = XNS(IN-1,ISURF)+XNS(IN,ISURF) TEMP2 = XNS(IN+1,ISURF)+XNS(IN,ISURF) TEMP3 = XNS(IN-1,ISURF)+XNS(IN,ISURF)*2.D0+XNS(IN+1,ISURF) TEMP4 = DS(IN-1,ISURF)+DS(IN,ISURF) VPORE(IN,ISURF) = 4.D0*(TEMP1*FLUX(IN-1,ISURF)- + TEMP2*FLUX(IN,ISURF)) + /(TEMP3*TEMP4) FLUX(1,2) = DDs*(AK(2,2)-AK(1,2))/DS(1,2) !!!opposite sign IN = 2 ISURF = 2 TEMP1 = XNS(IN-1,ISURF)+XNS(IN,ISURF) TEMP2 = XNS(IN+1,ISURF)+XNS(IN,ISURF) TEMP3 = XNS(IN-1,ISURF)+XNS(IN,ISURF)*2.D0+XNS(IN+1,ISURF) TEMP4 = DS(IN-1,ISURF)+DS(IN,ISURF) VPORE(IN,ISURF) = 4.D0*(TEMP1*FLUX(IN-1,ISURF)- + TEMP2*FLUX(IN,ISURF)) + /(TEMP3*TEMP4) RETURN END SUBROUTINE FILENAME(IFILE,JNAME,KJOB,DFILE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*11 JNAME CHARACTER*20 DFILE IF(IFILE.EQ.0) DFILE=JNAME(:KJOB)//'0.up' IF(IFILE.EQ.1) DFILE=JNAME(:KJOB)//'1.up' IF(IFILE.EQ.2) DFILE=JNAME(:KJOB)//'2.up' IF(IFILE.EQ.3) DFILE=JNAME(:KJOB)//'3.up' IF(IFILE.EQ.4) DFILE=JNAME(:KJOB)//'4.up' IF(IFILE.EQ.5) DFILE=JNAME(:KJOB)//'5.up' IF(IFILE.EQ.6) DFILE=JNAME(:KJOB)//'6.up' IF(IFILE.EQ.7) DFILE=JNAME(:KJOB)//'7.up' IF(IFILE.EQ.8) DFILE=JNAME(:KJOB)//'8.up' IF(IFILE.EQ.9) DFILE=JNAME(:KJOB)//'9.up' IF(IFILE.EQ.10) DFILE=JNAME(:KJOB)//'10.up' IF(IFILE.EQ.11) DFILE=JNAME(:KJOB)//'11.up' IF(IFILE.EQ.12) DFILE=JNAME(:KJOB)//'12.up' IF(IFILE.EQ.13) DFILE=JNAME(:KJOB)//'13.up' IF(IFILE.EQ.14) DFILE=JNAME(:KJOB)//'14.up' IF(IFILE.EQ.15) DFILE=JNAME(:KJOB)//'15.up' IF(IFILE.EQ.16) DFILE=JNAME(:KJOB)//'16.up' IF(IFILE.EQ.17) DFILE=JNAME(:KJOB)//'17.up' IF(IFILE.EQ.18) DFILE=JNAME(:KJOB)//'18.up' IF(IFILE.EQ.19) DFILE=JNAME(:KJOB)//'19.up' IF(IFILE.EQ.20) DFILE=JNAME(:KJOB)//'20.up' IF(IFILE.EQ.21) DFILE=JNAME(:KJOB)//'21.up' IF(IFILE.EQ.22) DFILE=JNAME(:KJOB)//'22.up' IF(IFILE.EQ.23) DFILE=JNAME(:KJOB)//'23.up' IF(IFILE.EQ.24) DFILE=JNAME(:KJOB)//'24.up' IF(IFILE.EQ.25) DFILE=JNAME(:KJOB)//'25.up' IF(IFILE.EQ.26) DFILE=JNAME(:KJOB)//'26.up' IF(IFILE.EQ.27) DFILE=JNAME(:KJOB)//'27.up' IF(IFILE.EQ.28) DFILE=JNAME(:KJOB)//'28.up' IF(IFILE.EQ.29) DFILE=JNAME(:KJOB)//'29.up' IF(IFILE.EQ.30) DFILE=JNAME(:KJOB)//'30.up' IF(IFILE.EQ.31) DFILE=JNAME(:KJOB)//'31.up' IF(IFILE.EQ.32) DFILE=JNAME(:KJOB)//'32.up' IF(IFILE.EQ.33) DFILE=JNAME(:KJOB)//'33.up' IF(IFILE.EQ.34) DFILE=JNAME(:KJOB)//'34.up' IF(IFILE.EQ.35) DFILE=JNAME(:KJOB)//'35.up' IF(IFILE.EQ.36) DFILE=JNAME(:KJOB)//'36.up' IF(IFILE.EQ.37) DFILE=JNAME(:KJOB)//'37.up' IF(IFILE.EQ.38) DFILE=JNAME(:KJOB)//'38.up' IF(IFILE.EQ.39) DFILE=JNAME(:KJOB)//'39.up' IF(IFILE.EQ.40) DFILE=JNAME(:KJOB)//'40.up' IF(IFILE.EQ.41) DFILE=JNAME(:KJOB)//'41.up' IF(IFILE.EQ.42) DFILE=JNAME(:KJOB)//'42.up' IF(IFILE.EQ.43) DFILE=JNAME(:KJOB)//'43.up' IF(IFILE.EQ.44) DFILE=JNAME(:KJOB)//'44.up' IF(IFILE.EQ.45) DFILE=JNAME(:KJOB)//'45.up' IF(IFILE.EQ.46) DFILE=JNAME(:KJOB)//'46.up' IF(IFILE.EQ.47) DFILE=JNAME(:KJOB)//'47.up' IF(IFILE.EQ.48) DFILE=JNAME(:KJOB)//'48.up' IF(IFILE.EQ.49) DFILE=JNAME(:KJOB)//'49.up' IF(IFILE.EQ.50) DFILE=JNAME(:KJOB)//'50.up' IF(IFILE.EQ.51) DFILE=JNAME(:KJOB)//'51.up' IF(IFILE.EQ.52) DFILE=JNAME(:KJOB)//'52.up' IF(IFILE.EQ.53) DFILE=JNAME(:KJOB)//'53.up' IF(IFILE.EQ.54) DFILE=JNAME(:KJOB)//'54.up' IF(IFILE.EQ.55) DFILE=JNAME(:KJOB)//'55.up' IF(IFILE.EQ.56) DFILE=JNAME(:KJOB)//'56.up' IF(IFILE.EQ.57) DFILE=JNAME(:KJOB)//'57.up' IF(IFILE.EQ.58) DFILE=JNAME(:KJOB)//'58.up' IF(IFILE.EQ.59) DFILE=JNAME(:KJOB)//'59.up' IF(IFILE.EQ.60) DFILE=JNAME(:KJOB)//'60.up' IF(IFILE.EQ.61) DFILE=JNAME(:KJOB)//'61.up' IF(IFILE.EQ.62) DFILE=JNAME(:KJOB)//'62.up' IF(IFILE.EQ.63) DFILE=JNAME(:KJOB)//'63.up' IF(IFILE.EQ.64) DFILE=JNAME(:KJOB)//'64.up' IF(IFILE.EQ.65) DFILE=JNAME(:KJOB)//'65.up' IF(IFILE.EQ.66) DFILE=JNAME(:KJOB)//'66.up' IF(IFILE.EQ.67) DFILE=JNAME(:KJOB)//'67.up' IF(IFILE.EQ.68) DFILE=JNAME(:KJOB)//'68.up' IF(IFILE.EQ.69) DFILE=JNAME(:KJOB)//'69.up' IF(IFILE.EQ.70) DFILE=JNAME(:KJOB)//'70.up' IF(IFILE.EQ.71) DFILE=JNAME(:KJOB)//'71.up' IF(IFILE.EQ.72) DFILE=JNAME(:KJOB)//'72.up' IF(IFILE.EQ.73) DFILE=JNAME(:KJOB)//'73.up' IF(IFILE.EQ.74) DFILE=JNAME(:KJOB)//'74.up' IF(IFILE.EQ.75) DFILE=JNAME(:KJOB)//'75.up' IF(IFILE.EQ.76) DFILE=JNAME(:KJOB)//'76.up' IF(IFILE.EQ.77) DFILE=JNAME(:KJOB)//'77.up' IF(IFILE.EQ.78) DFILE=JNAME(:KJOB)//'78.up' IF(IFILE.EQ.79) DFILE=JNAME(:KJOB)//'79.up' IF(IFILE.EQ.80) DFILE=JNAME(:KJOB)//'80.up' IF(IFILE.EQ.81) DFILE=JNAME(:KJOB)//'81.up' IF(IFILE.EQ.82) DFILE=JNAME(:KJOB)//'82.up' IF(IFILE.EQ.83) DFILE=JNAME(:KJOB)//'83.up' IF(IFILE.EQ.84) DFILE=JNAME(:KJOB)//'84.up' IF(IFILE.EQ.85) DFILE=JNAME(:KJOB)//'85.up' IF(IFILE.EQ.86) DFILE=JNAME(:KJOB)//'86.up' IF(IFILE.EQ.87) DFILE=JNAME(:KJOB)//'87.up' IF(IFILE.EQ.88) DFILE=JNAME(:KJOB)//'88.up' IF(IFILE.EQ.89) DFILE=JNAME(:KJOB)//'89.up' IF(IFILE.EQ.90) DFILE=JNAME(:KJOB)//'90.up' IF(IFILE.EQ.91) DFILE=JNAME(:KJOB)//'91.up' IF(IFILE.EQ.92) DFILE=JNAME(:KJOB)//'92.up' IF(IFILE.EQ.93) DFILE=JNAME(:KJOB)//'93.up' IF(IFILE.EQ.94) DFILE=JNAME(:KJOB)//'94.up' IF(IFILE.EQ.95) DFILE=JNAME(:KJOB)//'95.up' IF(IFILE.EQ.96) DFILE=JNAME(:KJOB)//'96.up' IF(IFILE.EQ.97) DFILE=JNAME(:KJOB)//'97.up' IF(IFILE.EQ.98) DFILE=JNAME(:KJOB)//'98.up' IF(IFILE.EQ.99) DFILE=JNAME(:KJOB)//'99.up' IF(IFILE.EQ.100) DFILE=JNAME(:KJOB)//'100.up' IF(IFILE.EQ.101) DFILE=JNAME(:KJOB)//'101.up' IF(IFILE.EQ.102) DFILE=JNAME(:KJOB)//'102.up' IF(IFILE.EQ.103) DFILE=JNAME(:KJOB)//'103.up' IF(IFILE.EQ.104) DFILE=JNAME(:KJOB)//'104.up' IF(IFILE.EQ.105) DFILE=JNAME(:KJOB)//'105.up' IF(IFILE.EQ.106) DFILE=JNAME(:KJOB)//'106.up' IF(IFILE.EQ.107) DFILE=JNAME(:KJOB)//'107.up' IF(IFILE.EQ.108) DFILE=JNAME(:KJOB)//'108.up' IF(IFILE.EQ.109) DFILE=JNAME(:KJOB)//'109.up' IF(IFILE.EQ.110) DFILE=JNAME(:KJOB)//'110.up' IF(IFILE.EQ.111) DFILE=JNAME(:KJOB)//'111.up' IF(IFILE.EQ.112) DFILE=JNAME(:KJOB)//'112.up' IF(IFILE.EQ.113) DFILE=JNAME(:KJOB)//'113.up' IF(IFILE.EQ.114) DFILE=JNAME(:KJOB)//'114.up' IF(IFILE.EQ.115) DFILE=JNAME(:KJOB)//'115.up' IF(IFILE.EQ.116) DFILE=JNAME(:KJOB)//'116.up' IF(IFILE.EQ.117) DFILE=JNAME(:KJOB)//'117.up' IF(IFILE.EQ.118) DFILE=JNAME(:KJOB)//'118.up' IF(IFILE.EQ.119) DFILE=JNAME(:KJOB)//'119.up' IF(IFILE.EQ.120) DFILE=JNAME(:KJOB)//'120.up' IF(IFILE.EQ.121) DFILE=JNAME(:KJOB)//'121.up' IF(IFILE.EQ.122) DFILE=JNAME(:KJOB)//'122.up' IF(IFILE.EQ.123) DFILE=JNAME(:KJOB)//'123.up' IF(IFILE.EQ.124) DFILE=JNAME(:KJOB)//'124.up' IF(IFILE.EQ.125) DFILE=JNAME(:KJOB)//'125.up' IF(IFILE.EQ.126) DFILE=JNAME(:KJOB)//'126.up' IF(IFILE.EQ.127) DFILE=JNAME(:KJOB)//'127.up' IF(IFILE.EQ.128) DFILE=JNAME(:KJOB)//'128.up' IF(IFILE.EQ.129) DFILE=JNAME(:KJOB)//'129.up' IF(IFILE.EQ.130) DFILE=JNAME(:KJOB)//'130.up' IF(IFILE.EQ.131) DFILE=JNAME(:KJOB)//'131.up' IF(IFILE.EQ.132) DFILE=JNAME(:KJOB)//'132.up' IF(IFILE.EQ.133) DFILE=JNAME(:KJOB)//'133.up' IF(IFILE.EQ.134) DFILE=JNAME(:KJOB)//'134.up' IF(IFILE.EQ.135) DFILE=JNAME(:KJOB)//'135.up' IF(IFILE.EQ.136) DFILE=JNAME(:KJOB)//'136.up' IF(IFILE.EQ.137) DFILE=JNAME(:KJOB)//'137.up' IF(IFILE.EQ.138) DFILE=JNAME(:KJOB)//'138.up' IF(IFILE.EQ.139) DFILE=JNAME(:KJOB)//'139.up' IF(IFILE.EQ.140) DFILE=JNAME(:KJOB)//'140.up' IF(IFILE.EQ.141) DFILE=JNAME(:KJOB)//'141.up' IF(IFILE.EQ.142) DFILE=JNAME(:KJOB)//'142.up' IF(IFILE.EQ.143) DFILE=JNAME(:KJOB)//'143.up' IF(IFILE.EQ.144) DFILE=JNAME(:KJOB)//'144.up' IF(IFILE.EQ.145) DFILE=JNAME(:KJOB)//'145.up' IF(IFILE.EQ.146) DFILE=JNAME(:KJOB)//'146.up' IF(IFILE.EQ.147) DFILE=JNAME(:KJOB)//'147.up' IF(IFILE.EQ.148) DFILE=JNAME(:KJOB)//'148.up' IF(IFILE.EQ.149) DFILE=JNAME(:KJOB)//'149.up' IF(IFILE.EQ.150) DFILE=JNAME(:KJOB)//'150.up' IF(IFILE.EQ.151) DFILE=JNAME(:KJOB)//'151.up' IF(IFILE.EQ.152) DFILE=JNAME(:KJOB)//'152.up' IF(IFILE.EQ.153) DFILE=JNAME(:KJOB)//'153.up' IF(IFILE.EQ.154) DFILE=JNAME(:KJOB)//'154.up' IF(IFILE.EQ.155) DFILE=JNAME(:KJOB)//'155.up' IF(IFILE.EQ.156) DFILE=JNAME(:KJOB)//'156.up' IF(IFILE.EQ.157) DFILE=JNAME(:KJOB)//'157.up' IF(IFILE.EQ.158) DFILE=JNAME(:KJOB)//'158.up' IF(IFILE.EQ.159) DFILE=JNAME(:KJOB)//'159.up' IF(IFILE.EQ.160) DFILE=JNAME(:KJOB)//'160.up' IF(IFILE.EQ.161) DFILE=JNAME(:KJOB)//'161.up' IF(IFILE.EQ.162) DFILE=JNAME(:KJOB)//'162.up' IF(IFILE.EQ.163) DFILE=JNAME(:KJOB)//'163.up' IF(IFILE.EQ.164) DFILE=JNAME(:KJOB)//'164.up' IF(IFILE.EQ.165) DFILE=JNAME(:KJOB)//'165.up' IF(IFILE.EQ.166) DFILE=JNAME(:KJOB)//'166.up' IF(IFILE.EQ.167) DFILE=JNAME(:KJOB)//'167.up' IF(IFILE.EQ.168) DFILE=JNAME(:KJOB)//'168.up' IF(IFILE.EQ.169) DFILE=JNAME(:KJOB)//'169.up' IF(IFILE.EQ.170) DFILE=JNAME(:KJOB)//'170.up' IF(IFILE.EQ.171) DFILE=JNAME(:KJOB)//'171.up' IF(IFILE.EQ.172) DFILE=JNAME(:KJOB)//'172.up' IF(IFILE.EQ.173) DFILE=JNAME(:KJOB)//'173.up' IF(IFILE.EQ.174) DFILE=JNAME(:KJOB)//'174.up' IF(IFILE.EQ.175) DFILE=JNAME(:KJOB)//'175.up' IF(IFILE.EQ.176) DFILE=JNAME(:KJOB)//'176.up' IF(IFILE.EQ.177) DFILE=JNAME(:KJOB)//'177.up' IF(IFILE.EQ.178) DFILE=JNAME(:KJOB)//'178.up' IF(IFILE.EQ.179) DFILE=JNAME(:KJOB)//'179.up' IF(IFILE.EQ.180) DFILE=JNAME(:KJOB)//'180.up' IF(IFILE.EQ.181) DFILE=JNAME(:KJOB)//'181.up' IF(IFILE.EQ.182) DFILE=JNAME(:KJOB)//'182.up' IF(IFILE.EQ.183) DFILE=JNAME(:KJOB)//'183.up' IF(IFILE.EQ.184) DFILE=JNAME(:KJOB)//'184.up' IF(IFILE.EQ.185) DFILE=JNAME(:KJOB)//'185.up' IF(IFILE.EQ.186) DFILE=JNAME(:KJOB)//'186.up' IF(IFILE.EQ.187) DFILE=JNAME(:KJOB)//'187.up' IF(IFILE.EQ.188) DFILE=JNAME(:KJOB)//'188.up' IF(IFILE.EQ.189) DFILE=JNAME(:KJOB)//'189.up' IF(IFILE.EQ.190) DFILE=JNAME(:KJOB)//'190.up' IF(IFILE.EQ.191) DFILE=JNAME(:KJOB)//'191.up' IF(IFILE.EQ.192) DFILE=JNAME(:KJOB)//'192.up' IF(IFILE.EQ.193) DFILE=JNAME(:KJOB)//'193.up' IF(IFILE.EQ.194) DFILE=JNAME(:KJOB)//'194.up' IF(IFILE.EQ.195) DFILE=JNAME(:KJOB)//'195.up' IF(IFILE.EQ.196) DFILE=JNAME(:KJOB)//'196.up' IF(IFILE.EQ.197) DFILE=JNAME(:KJOB)//'197.up' IF(IFILE.EQ.198) DFILE=JNAME(:KJOB)//'198.up' IF(IFILE.EQ.199) DFILE=JNAME(:KJOB)//'199.up' IF(IFILE.EQ.200) DFILE=JNAME(:KJOB)//'200.up' IF(IFILE.GT.200) DFILE=JNAME(:KJOB)//'_last.up' RETURN END SUBROUTINE CLEAR(A,N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(N) DO 100 I=1,N A(I)=0.0D0 100 CONTINUE RETURN END CENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND