C**************************************************************** C CAMARAS DE MAGMA RODEADA DE ROCA C MODELO NO NEWTONIANOS PARA MAGMA EN CONVECCIÓN NATURAL C**************************************************************** BLOCK DATA CHARACTER*10 TITLE LOGICAL LSOLVE,LPRINT,LBLK,LSTOP COMMON/CHAR/TITLE(13) COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, 1 IST,JST,ITER,LAST,RELAX(13),TTIME,DT,XL,YL, 2 IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON COMMON/CNTL/LSTOP COMMON/SORC/SMAX,SSUM COMMON/COEF/FLOW,DIFF,ACOF DATA NFMAX,NP,NRHO,NGAM/10,11,12,13/ DATA LSTOP/1*.FALSE./ DATA LBLK/10*.TRUE./ DATA MODE,TTIME,ITER/1,0.,0/ DATA RELAX(13)/1./ DATA NTIMES/10*1/ DATA DT,IPREF,JPREF,RHOCON/1.E+10,1,1,1./ C**************************************************************** C ENTRADA DE DATOS A ETIQUETA COMUN USANDO DATOS DECLARADOS C *********************************************************** DATA TITLE(1),TITLE(2),TITLE(3),TITLE(4),TITLE(5),TITLE(6), %TITLE(7),TITLE(11)/ 1'VEL U',' VEL V',' STR FN',' TEMP','SOL','GAMAP','ETA','PRESSURE'/ DATA RELAX(1),RELAX(2),RELAX(4),RELAX(5),RELAX(6),RELAX(11) 1/0.25,0.25,1.,1.,1.,0.8/ DATA (LSOLVE(I),LPRINT(I),I=1,5),LPRINT(6),LPRINT(7), %LPRINT(11)/13*.TRUE./ DATA LAST/1e9/ END C**************************************************************** C DEFINICIÓN DE LOS ARCHIVOS DE DATOS C**************************************************************** SUBROUTINE OPENF OPEN(UNIT=3,FILE='prueba.plt') OPEN(UNIT=4,FILE='COEF_NO_NEWT.TXT') OPEN(UNIT=1504,FILE='Todo_TpoFinal.PLT') OPEN(UNIT=3,FILE='prueba.plt') OPEN(UNIT=1505,FILE='TodoAdmim.PLT') OPEN(UNIT=1506,FILE='TodoDimensional.PLT') OPEN(UNIT=1507,FILE='PerfilX_U.PLT') OPEN(UNIT=1607,FILE='PerfilX_V.PLT') OPEN(UNIT=1508,FILE='PerfilX_Temp.PLT') OPEN(UNIT=1608,FILE='PerfilX.ConcPLT') OPEN(UNIT=1509,FILE='PerfilY_U.PLT') OPEN(UNIT=1609,FILE='PerfilY_V.PLT') OPEN(UNIT=1510,FILE='PerfilY_Temp.PLT') OPEN(UNIT=1610,FILE='PerfilY_Conc.PLT') OPEN(UNIT=70,FILE='ConcDim.plt') OPEN(UNIT=71,FILE='TempDim.plt') OPEN(UNIT=72,FILE='VelocDim.plt') OPEN(UNIT=73,FILE='TPROM.PLT') OPEN(UNIT=74,FILE='CPROM.PLT') OPEN(UNIT=75,FILE='UPROM.PLT') OPEN(UNIT=76,FILE='VPROM.PLT') OPEN(UNIT=77,FILE='etaPROM.PLT') OPEN(UNIT=78,FILE='eta.PLT') OPEN(UNIT=88,FILE='VresProm.PLT') OPEN(UNIT=101,FILE='DURACION.TXT') OPEN(UNIT=2402,FILE='PUNTO_MAS_CALIENTE.DAT') RETURN END C**************************************************************** SUBROUTINE USER LOGICAL LSOLVE,LPRINT,LBLK,LSTOP PARAMETER (NOD=500) !NUMERO MAXIMO DE NODOS COMMON F(NOD,NOD,10),P(NOD,NOD),RHO(NOD,NOD),GAM(NOD,NOD), 1CON(NOD,NOD), 1 AIP(NOD,NOD),AIM(NOD,NOD),AJP(NOD,NOD),AJM(NOD,NOD),AP(NOD,NOD), 2 X(NOD),XU(NOD),XDIF(NOD),XCV(NOD),XCVS(NOD), 3 Y(NOD),YV(NOD),YDIF(NOD),YCV(NOD),YCVS(NOD), 4 YCVR(NOD),YCVRS(NOD),ARX(NOD),ARXJ(NOD),ARXJP(NOD), 5 R(NOD),RMN(NOD),SX(NOD),SXMN(NOD),XCVI(NOD),XCVIP(NOD) COMMON DU(NOD,NOD),DV(NOD,NOD),FV(NOD),FVP(NOD), 1 FX(NOD),FXM(NOD),FY(NOD),FYM(NOD),PT(NOD),QT(NOD) COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, 1 IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, 2 IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCO COMMON/CNTL/LSTOP COMMON/SORC/SMAX,SSUM COMMON/COEF/FLOW,DIFF,ACOF DIMENSION U(NOD,NOD),V(NOD,NOD),ABC1(NOD,NOD),ABC2(NOD,NOD), 1EFG1(NOD,NOD),EFG2(NOD,NOD),CONVER1(NOD,NOD),CONVER2(NOD,NOD), 2GHI1(NOD,NOD),GHI2(NOD,NOD),CONVER3(NOD,NOD), 3JKL1(NOD,NOD),JKL2(NOD,NOD),CONVER4(NOD,NOD), 4UANT(NOD,NOD),VANT(NOD,NOD),TANT(NOD,NOD),CANT(NOD,NOD) EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),STR(1,1)) DIMENSION TH(NOD),THU(NOD),THDIF(NOD),THCV(NOD),THCVS(NOD) EQUIVALENCE(X,TH),(XU,THU),(XDIF,THDIF),(XCV,THCV), 1(XCVS,THCVS),(XL,THL) COMMON/MVAT/NTYV,NTXV,VTYV(30),VTXV(30),YLTV(30),XLTV(30),ITRAM, 1COORD,RX(30),RY(30) C**************************************************************** DIMENSION T(NOD,NOD),C(NOD,NOD),STR(NOD,NOD),RA(NOD,NOD), 1PR(NOD,NOD),dE(NOD,NOD) EQUIVALENCE(F(1,1,4),T(1,1)),(F(1,1,5),C(1,1)) real ETAI(NOD,NOD) REAL TA(2) real cs1dim(nod),cs2dim(nod),cs3dim(nod),cs4dim(nod) &,cs5dim(nod),cs6dim(nod) &,ci3dim(nod),ci4dim(nod),ci5dim(nod),ci6dim(nod),ci7dim(nod) REAL TFILM,etain,varaux1 &,varaux2,varaux3,RA0,PR0 real alfaStar,etaprom real FuenteVel,fuenteT,fuenteC,fuenteTS,fuenteCS &,fuenteVelFlota,GamModelV,GamModelTS,GamModelT,GamModelC,le &,deltaTad REAL xlong,ylong,tdimen,lo,xxl,yyl,LL,TempAdI,TempAdF,LongC,LongCL integer xmin1,xmin2,xmin3,xmin4,xmin5,xmax1,xmax2,xmax3,xmax4 &,ymax1,ymax2,ymax3,ymax4,ymax5,ymax6,ymax7,ymax8,ymax9,ymax10 INTEGER CONTI,TINIC,TEX,xnodo,ynodo,salida,imprimir1 &,yprof1,xprof1,imprime integer PARTES(NOD,NOD) real camara(nod,nod) C**************************************************************** C**************************************************************** ENTRY GRID MODE=1 ! goto 123 xxl = 24000 LL = 24000 !longitud mayor en metros yyl = 6000 LongC = 3000 !altura de la camara NTXV = 5 NTYV = 5 C-------CANTIDAD DE NODOS EN X--------(198 NODOS) VTXV(1)=40 VTXV(2)=20 VTXV(3)=76 VTXV(4)=20 !22 VTXV(5)=40 !24 ! VTXV(6)=7 ! VTXV(7)=56 xnodo = vtxV(1) + vtxV(2) + vtxV(3) + vtxV(4) + vtxV(5) + 2 write(*,*)'nodos en zona en x', xnodo C-------CANTIDAD DE NODOS EN Y--------(60 NODOS) VTYV(1)=15 VTYV(2)=15 VTYV(3)=16 VTYV(4)=5 VTYV(5)=5 ynodo = vtyV(1) + vtyV(2) + vtyV(3) + vtyV(4) + vtyV(5) + 2 write(*,*)'nodos en zona en y', ynodo C-------LONGITUD DE LA ZONA EN X------ XLTV(1)=6800/LongC !2.8 XLTV(2)=1000/LongC !0.2 XLTV(3)=9000/LongC !0.685714 XLTV(4)=1000/LongC !0.628571 XLTV(5)=6200/LongC !0.685714 ! XLTV(6)=0.2 ! XLTV(7)=2.8 xlong = XLTV(1) + XLTV(2) + XLTV(3) + XLTV(4) + XLTV(5) write(*,*)'longitud zona en x', xlong C-------LONGITUD DE LA ZONA EN Y------ YLTV(1)=1000/LongC YLTV(2)=1400/LongC YLTV(3)=1000/LongC YLTV(4)=1350/LongC YLTV(5)=1250/LongC ylong = yLTV(1) + yLTV(2) + yLTV(3) + yLTV(4) + yLTV(5) write(*,*)'longitud zona en y', ylong*longc C-------REFINACION EN LA DIRECCION X-- RX(1)=1. RX(2)=1. RX(3)=1. RX(4)=1. RX(5)=1. ! RX(6)=1. ! RX(7)=1. C-------REFINACION EN LA DIRECCION Y-- RY(1)=1. RY(2)=1. RY(3)=1. RY(4)=1. RY(5)=1. CALL MALLA_VAT 123 continue ! xxl = 24000 ! LL = 24000 !longitud caracteristica en metros ! yyl = 6000 ! LongC = 3000 !longitud característica en metros ! L1 = 198 ! XL = xxl/LongC ! M1 = 60 ! YL = yyl/LongC ! call ugrid RETURN C**************************************************************** C**************************************************************** ENTRY START DTPO = 1.e1 TIEMPO = 0. TFIN = 1e8 !TIEMPO FINAL pero para si alcanza temperaturas promedio bajo los 1448 TIEMPOFINAL=TIMEF() !TIEMPO COMPUTACIONAL EMAX=0. ETMAX=0. ITERMIN=7.;ITERMAX=30. !ITERACIÓN MÍNIMA Y MAXIMA CRITERIO=1.E-6 !criterio de convergencia para paso de tiempo TPRINT=0. TOUT=1e-2 !*************************************************************** salida = 0 imprime = 0 SMALL=1.E-30 partes = 1 !roca anfitriona xmin1=0 xmax1=0 xmax2=0 ymin1=0 ymin2=0 ymax1=0 ymax2=0 !forma rectangular ! goto 101 ! do i = 1,l1 ! if(x(i).le.((xxl/2.)-3000.)/LongC)xmin1 = i ! if(x(i).le.((xxl/2.)+3000.)/LongC)xmax1 = i ! enddo ! do j = 1,m1 ! if(y(j).le.(yyl-3000.)/LongC) ymax1 = j ! enddo ! DO 101 J=1,M1 ! DO 101 I=1,L1 ! if(i.ge.xmin1.and.i.le.xmax1.and.j.le.ymax1)then ! partes(i,j)=0 ! endif ! 101 CONTINUE !******forma eliptica do j = 1,m1 if(y(j).le.(2211)/LongC) ymax1 = j if(y(j).le.(2766)/LongC) ymax2 = j if(y(j).le.(3000)/LongC) ymax3 = j if(y(j).le.(2935)/LongC) ymax4 = j if(y(j).le.(2768)/LongC) ymax5 = j if(y(j).le.(2769)/LongC) ymax6 = j if(y(j).le.(2627)/LongC) ymax7 = j if(y(j).le.(2661)/LongC) ymax8 = j if(y(j).le.(2318)/LongC) ymax9 = j if(y(j).le.(2483)/LongC) ymax10 = j if(y(j).le.(1.)/LongC) ymin1 = j if(y(j).le.(2100)/LongC) ymin2 = j if(y(j).le.(521)/LongC) ymin3 = j if(y(j).le.(796)/LongC) ymin4 = j if(y(j).le.(1073)/LongC) ymin5 = j if(y(j).le.(1312)/LongC) ymin6 = j if(y(j).le.(1516)/LongC) ymin7 = j if(y(j).le.(1700)/LongC) ymin8 = j if(y(j).le.(1960)/LongC) ymin9 = j enddo do i = 1,l1 if(x(i).le.(7000)/LongC)xmin1 = i if(x(i).le.(12257)/LongC)xmin2 = i if(x(i).le.(13885)/LongC)xmin3 = i if(x(i).le.(15428)/LongC)xmin4 = i if(x(i).le.(16826)/LongC)xmin5 = i if(x(i).le.(7000)/LongC)xmax1 = i if(x(i).le.(11795)/LongC)xmax2 = i if(x(i).le.(11794)/LongC)xmax3 = i if(x(i).le.(13133)/LongC)xmax4 = i if(x(i).le.(13389)/LongC)xmax5 = i if(x(i).le.(14460)/LongC)xmax6 = i if(x(i).le.(15836)/LongC)xmax7 = i if(x(i).le.(18000)/LongC)xmax8 = i enddo !******measure points for the profiles do j = 1,m1 if(y(j).le.(yyl-4500)/LongC) yprof1 = j enddo do i = 1,l1 if(x(i).le.(xxl/2.)/LongC) xprof1 = i enddo !***** arriba izquierda cs1dim = 0. do j = ymax1,ymax2 do i = xmax1,xmax2 cs1dim(i) = (-2.79104210E-18*(x(i)*LongC)**6 & + 1.60078795E-13*(x(i)*LongC)**5 & - 3.80443891E-09*(x(i)*LongC)**4 & + 4.79537286E-05*(x(i)*LongC)**3 & - 3.38112514E-01*(x(i)*LongC)**2 & + 1.26459526E+03*(x(i)*LongC) - 1.95815530E+06)/LongC if(y(j).le.cs1dim(i))then partes(i,j) = 0 endif enddo enddo !***** arriba medio cs6dim = 0. do j = ymin8,ymax3 do i = xmax2,xmax4 cs6dim(i) = 3000/LongC if(y(j).le.cs6dim(i))then partes(i,j) = 0 endif enddo enddo !***** arriba derecha 1 cs2dim = 0. do j = ymax9,ymax10 do i = xmax7,xmax8-4 cs2dim(i) = (-6.28785583E-14*(x(i)*LongC)**5 & + 5.31719532E-09*(x(i)*LongC)**4 & - 1.79757534E-04*(x(i)*LongC)**3 & + 3.03683503E+00*(x(i)*LongC)**2 & - 2.56378386E+04*(x(i)*LongC) + 8.65302167E+07)/LongC if(y(j).le.cs2dim(i))then partes(i,j) = 0 endif enddo enddo !***** arriba derecha 2 cs3dim = 0. do j = ymin8,ymax8 do i = xmax6,xmax7 cs3dim(i) = (-7.97247969E-08*(x(i)*LongC)**3 & + 3.64603962E-03*(x(i)*LongC)**2 & - 5.55663080E+01*(x(i)*LongC) + 2.84836825E+05)/LongC if(y(j).le.cs3dim(i))then partes(i,j) = 0 endif enddo enddo !***** arriba derecha 3 cs4dim = 0. do j = ymin8,ymax6 do i = xmax5,xmax6 cs4dim(i) = (-2.92472483E-16*(x(i)*LongC)**5 & + 2.05562047E-11*(x(i)*LongC)**4 & - 5.77715953E-07*(x(i)*LongC)**3 & + 8.11538958E-03*(x(i)*LongC)**2 & - 5.69796568E+01*(x(i)*LongC) + 1.62731539E+05)/LongC if(y(j).le.cs4dim(i))then partes(i,j) = 0 endif enddo enddo !***** arriba derecha 4 cs5dim = 0. do j = ymin8,ymax4 do i = xmax4,xmax5 cs5dim(i) = (2935)/LongC if(y(j).le.cs5dim(i))then partes(i,j) = 0 endif enddo enddo !**** abajo izquierda 1 ci3dim = 0. do j = ymin1,ymin2 do i = xmin1,xmin2 ci3dim(i)= (-2.38353227E-15*(x(i)*LongC)**5 & + 1.15020924E-10*(x(i)*LongC)**4 & - 2.21004236E-06*(x(i)*LongC)**3 & + 2.11177530E-02*(x(i)*LongC)**2 & - 1.00556013E+02*(x(i)*LongC) + 1.92867532E+05)/LongC if(y(j).ge.ci3dim(i))then partes(i,j) = 0 endif enddo enddo !**** abajo derecha 2 ci4dim = 0. do j = ymin1,ymin9 do i = xmin2,xmin3 ci4dim(i)= (1.46942793E-15*(x(i)*LongC)**5 & - 5.67835653E-11*(x(i)*LongC)**4 & + 4.20558880E-07*(x(i)*LongC)**3 & + 8.89814265E-03*(x(i)*LongC)**2 & - 1.55062531E+02*(x(i)*LongC) + 6.64491198E+05)/LongC if(y(j).ge.ci4dim(i))then partes(i,j) = 0 endif enddo enddo !**** abajo derecha 3 do j = ymin4,ymin9 do i = xmin3,xmin4 ci5dim(i)= (-3.53256839E-14*(x(i)*LongC)**5 & + 2.54566190E-09*(x(i)*LongC)**4 & - 7.33312702E-05*(x(i)*LongC)**3 & + 1.05550712E+00*(x(i)*LongC)**2 & - 7.59103964E+03*(x(i)*LongC) + 2.18218376E+07)/LongC if(y(j).ge.ci5dim(i))then partes(i,j) = 0 endif enddo enddo !**** abajo derecha 4 do j = ymin6,ymax1 do i = xmin4,xmin5 ci6dim(i)= (2.03590921E-14*(x(i)*LongC)**5 & - 1.56814913E-09*(x(i)*LongC)**4 & + 4.82086615E-05*(x(i)*LongC)**3 & - 7.39261408E-01*(x(i)*LongC)**2 & + 5.65354471E+03*(x(i)*LongC) - 1.72449199E+07)/LongC if(y(j).ge.ci6dim(i))then partes(i,j) = 0 endif enddo enddo ! goto 4321 !**** abajo derecha 5 do j = ymin8,ymax1 do i = xmin5,xmax8-4 ci7dim(i)= (1.49845254E-01*(x(i)*LongC) - 7.41759923E+02)/LongC ! write(*,*)x(i)*longc,ci7dim(i)*longc,y(j)*longc ! pause if(y(j).ge.ci7dim(i))then partes(i,j) = 0 endif enddo enddo 4321 continue camara = 1 volumen = 0. DO 401 J=1,M1 DO 401 I=1,L1 if(partes(i,j).eq.0)then volumen = volumen + xcv(i)*LongC*ycv(j)*LongC*1 camara(i,j) = 0 endif 401 CONTINUE write(*,*) 'volume chamber m2=',volumen,' km2=',volumen/1.e6 pause ! WRITE(3,*) 'TITLE = "prueba1"' ! WRITE(3,*) 'VARIABLES="X","Y","partes"' ! WRITE(3,*) 'ZONE T="',TIEMPO,'" I=',L1,', J=',M1 ! DO J=1,M1 ! DO I=1,L1 ! WRITE(3,*) X(I)*longc,Y(J)*longc,partes(I,J) ! enddo ! enddo ! stop !***********************LIMITES DE TEMPERAURA FLUIDO*********************************** TINIC = 1498 !TEMPERATURA INICIAL magma TEX = 1448 !TEMPERATURA FINAL magma TFILM = (TEX+TINIC)/2. !temperatura media magma !*********************LIMITES dimensionales DE TEMPERATURA HOST ROCK******************* TGrad = 25 !varia 25 °C cada kilómetro TO = 25 + 273 !temperatura inicial menor en host rock TD = TGrad * (yyl/1000) +25 +273 !final mayor en host rock !**************************Temperaturas adimensionales********************************* TempAdI = (Tinic - To)/ (Tinic-To) !temp inicial adimensional magma TempAdF = (Tex - To) / (Tinic-To) !temp final adimensional magma toad = (to - To) / (Tinic-To) !temp arriba adimensional host Tdad = (td - To) / (Tinic-To) !temp abajo adimensional host PEND = -(Tdad-Toad) / (yl-0) !pendiente adimensional Taf = (Tinic + Tex) / 2 TaH = (TO + TD) / 2 ! write(*,*) tempadi,(Tinic - To)*tempadi+to ! pause !***********************LIMITES DE CONCENTRACIÓN FLUIDO*********************************** CINIC = 0.7 !CONCENTRACIÓN INICIAL magma riolitico SI02 CEX = 0.47 !CONCENTRACIÓN INTRUSIÓN magma basaltico SIO2 CFILM = (CEX+CINIC)/2. !CONCENTRACIÓN media magma !*********************LIMITES dimensionales DE CONCENTRACION HOST ROCK******************* CGrad = 0.001 !varia XXXX cada kilómetro CO = 0.4 !CONCENTRACIÓN silica inicial menor en host rock CD = CGrad * (yyl/1000)+CO !CONCENTRACIÓN FINAL mayor en host rock !**************************CONCENTRACIÓN adimensionales********************************* CempAdI = (Cinic - Co)/ (Cinic-Co) !CONC inicial adimensional magma CempAdF = (Cex - Co) / (Cinic-Co) !CONC final adimensional magma Coad = (Co - Co) / (Cinic-Co) !CONC arriba adimensional host Cdad = (Cd - Co) / (Cinic-Co) !CONC abajo adimensional host CEND = -(Cdad-Coad)/ (yl-0) !pendiente adimensional ! write(*,*) td,to,CempAdi ! pause !**************************lARGO CARACTERISTICO**************************************** LO = 6000 !LONGITUD m R0 = 0.637 !Longitud caracteristica para rayleigh en metros CaRRIGAN paper GE = 9.81 !GRAVEDAD CTEGAS = 8.31451 !GASES J/mol K !**************CTE Y PENDIENTE DE LA ENERGIA DE ACTIVACION VISCOSIDAD****************** YY0 = 1492896.7965 AA0 = -994.4827 !********************* PROPIEDADES DE HOST ROCK**************************************** DENSH = 2670. !densidad kg/m3 CPH = 1000. !calor específico J/kg K KH = 2.65 !conductividad termica W/°C m2 ALFAH = 9.93E-7 !difusividad term m2/s !************************ PENDIENTE DE GRADIENTE DE TEMPERATURA para host rock***** !******************* FLUJO DE CALOR SOBRE Y POR DEBAJO****************************** Q1 = 0.065 QM1 = -0.065 !*******************parametros de adimensionalización en km****************************** KX = 3. !altura de camara MY = 3. !altura de camar en y KY = -6. !ancho de la cámara C**************** liquido interior ********************** DENS = 2600 !densidad magma basaltico CP = 1450 !calor específico magma COND = 0.6 !conductividad magma BETA = 5E-5 !coeficiente de expansión volumétrica etain = 100 !382.96 !viscosidad inicial de referencia magma basáltico Pa s ETAREF = ETAin !VISCOCIDAD DE REFERENCIA PA S eta0 = 15 !indice de consistencia K (Pa s) ALFA = COND/(DENS*CP) !difusividad termica fluido LongCL = 0.404 !espesor medio de la capa límite termica carrigan deltaT = 1.01e-2 !paper carrigan calculado en excel con datos de este artículo deltaTAd= (deltaT - To) / (Tinic-To) U0 = 8.9e-4 !MAXIMA VELOCIDAD EN EL FLUIDO paper Clark U0 = sqrt(ge*beta*deltaT*LongCL) !Maxima Veloc Fluido calculada PR0 = (ETAREF*cp)/cond !número de prandlt ! PR0 = 4E4 RA0 = (BETA*GE*deltaT*(LongC)**3) !NUMERO DE RAYLEIGH INICIAL en el fluido & /((ETAREF/dens)*ALFA) ! RA0 = 1e10 DIJ = 3E-8 !difusividad del Si en la magma m2/s LE = ALFA/DIJ !numero de lewis alfaStar = alfa/alfaH write(*,*)'Ra=',Ra0,' lewis=',Le,' Pr=',Pr0 write(*,*)'volumen m2=',Volumen,' km2=',volumen/1.e6 &,' velocidad maxima U0=',U0 ! &,cempadf,(cinic - co)*cempadf+co,coad pause !VALORES INICIALES******************************* DO 100 J=1,M1 DO 100 I=1,L1 U(I,J) = 0. V(I,J) = 0. t(i,j) = 0. c(i,j) = 0. !*******ENERGIA DE ACTIVACION dimensional solo magma****************** dE(I,J)=YY0+AA0*abs(((TINIC-TO)*T(I,J))+TO) !INICIALIZACIÓN DE RA Y PR *********************************** RA(I,J) = RA0 PR(I,J) = PR0 100 CONTINUE !*****condiciones iniciales en los bordes y dominio**************** DO 102 J=1,M1 DO 102 I=1,L1 IF(partes(I,J).EQ.0) AN=0.58925 !no newtoniano liquido interior 102 CONTINUE DO 302 I=1,L1 DO 302 J=1,M1 IF(partes(I,J).EQ.0)then T(I,J) = TempAdI !temp inicial en magma C(I,J) = CempAdi ! Si inicial en magma riolita else T(I,J)=TDad+PEND*(Y(J)) !temp inicial roca anfitriona C(I,J)= CDad+CEND*(Y(J)) !conc inicial roca ENDIF T(I,1)=TDad+PEND*(Y(1)) !temp inicial roca anfitriona U(I,J)=0. V(I,J)=0. 302 CONTINUE DO I=1,L1 DO J=1,M1 IF(partes(I,J).EQ.0.and.y(j)*LongC.le.1000)then partes(i,j) = 2 C(I,J) = CempAdF !Si basalto intruye en riolita AN = 0.58925 ENDIF enddo enddo ! WRITE(3,*) 'TITLE = "prueba1"' ! WRITE(3,*) 'VARIABLES="X","Y","TEMP","CONC"' ! WRITE(3,*) 'ZONE T="',TIEMPO,'" I=',L1,', J=',M1 ! DO J=1,M1 ! DO I=1,L1 ! WRITE(3,*) X(I)*LongC,Y(J)*LongC-6000, ! & partes(I,J),(cinic - co)*c(i,j)+co ! enddo ! enddo ! stop !**************************************************************** DO 110 J=1,M1 DO 110 I=1,L1 UANT(I,J)=U(I,J) VANT(I,J)=V(I,J) TANT(I,J)=T(I,J) CANT(I,J)=C(I,J) ABC1(I,J)=1. EFG1(I,J)=1. GHI1(I,J)=1. JKL1(I,J)=1 110 CONTINUE RETURN C**************************************************************** C**************************************************************** C************!solo lo hace para u y v velocidades**************** ENTRY DENSE RETURN C**************************************************************** C**************************************************************** C**************************************************************** ENTRY BOUND DO 300 I=1,L1 DO 300 J=1,M1 t(1,j) = t(2,j) t(l1,j) = t(l2,j) t(i,1) = t(i,2) c(i,1) = c(i,2) ! c(i,1) = c(i,2) IF(camara(I,J).EQ.1)THEN ! T(I,1)=T(I,2)+((YDIF(2)*Q1)/KH) C(I,J)=CDad+CEND*(Y(J)) !conc inicial roca ENDIF !VELOCIDADES CERO EN LAS PAREDES IF(partes(I,J).EQ.1)U(I,J)=0 IF(partes(I,J).EQ.1)V(I,J)=0 300 CONTINUE RETURN C**************************************************************** C**************************************************************** ENTRY OUTPUT CONT=0 DO 521 I=1,L1 DO 521 J=1,M1 ABC2(I,J)=U(I,J) EFG2(I,J)=V(I,J) GHI2(I,J)=T(I,J) JKL2(I,J)=C(I,J) CONVER1(I,J)=ABS(ABC2(I,J)-ABC1(I,J)) CONVER2(I,J)=ABS(EFG2(I,J)-EFG1(I,J)) CONVER3(I,J)=ABS(GHI2(I,J)-GHI1(I,J)) CONVER4(I,J)=ABS(JKL2(I,J)-JKL1(I,J)) ABC1(I,J)=ABC2(I,J) EFG1(I,J)=EFG2(I,J) GHI1(I,J)=GHI2(I,J) JKL1(I,J)=JKL2(I,J) !criterio de convergencia para impresión de resultados y proximo paso de tpo IF(CONVER1(I,J).LE.CRITERIO.AND.CONVER2(I,J).LE.CRITERIO.AND. 1CONVER3(I,J).LE.CRITERIO.AND.CONVER4(I,J).LE.CRITERIO)THEN CONT=CONT+1 ENDIF 521 CONTINUE 333 FORMAT(I6.2,F14.2,F14.2,2F10.4,5E12.4) WRITE(*,333)ITERN,TIEMPO,tdimen/3.154e7,tprom,cprom,vprom*U0 &,uprom*U0,((TINIC-TO)*t(xprof1,yprof1))+TO &,smax,VresProm*U0 !etai(xprof1,yprof1) !*U0 ! write(*,*) MaxTemp,(((TINIC-TO)*MaxTemp)+TO) itermin = 7 itermax = 33 IF(SMAX.Gt.1.E-3.and.dtpo.ge.1e-2 &.and.itern.ge.itermin) DTPO=DTPO*0.9 dtpo=10 ! if(smax.ge.1e-5)then ! relax(1)=0.075 ! relax(2)=0.055 ! relax(4)=0.8 ! relax(5)=0.8. ! relax(11)=0.8 ! else ! relax(1)=0.25 ! relax(2)=0.25 ! relax(4)=1. ! relax(5)=1. ! relax(11)=0.8 ! ! endif !INICIA IMPRESIÓN DE RESULTADOS Y PASO DE TIEMPO IF (CONT.GE.L1*M1*0.96.AND.ITERN.GE.ITERMIN.AND.SMAX.LE.1.E-4 !INTRODUCIR CRITERIO DE CONVERGENCIA %.OR.ITERN.EQ.itermax)THEN C****PARA CALCULAR variación maxima de la variable y aumentar paso de tiempo en los primeros segundos ! ETMAX=0. ! DO 406 I=1,L1 ! DO 406 J=1,M1 ! IF(ABS(U(I,J)-UANT(I,J)).GT.ETMAX)ETMAX=ABS(U(I,J)-UANT(I,J)) ! IF(ABS(V(I,J)-VANT(I,J)).GT.ETMAX)ETMAX=ABS(V(I,J)-VANT(I,J)) ! IF(ABS(T(I,J)-TANT(I,J)).GT.ETMAX)ETMAX=ABS(T(I,J)-TANT(I,J)) ! IF(ABS(C(I,J)-CANT(I,J)).GT.ETMAX)ETMAX=ABS(C(I,J)-CANT(I,J)) ! 406 CONTINUE IF(sMAX.Le.1E-4.AND.DTPO.LT.10)DTPO=DTPO*1.05 tdimen = (tiempo * LongC)/U0 varaux1 = 10 varaux2 = 100 varaux3 = 500 varaux4 = 1000 varaux5 = 2.5e3 varaux6 = 3e3 varaux7 = 3.5e3 varaux8 = 4e3 varaux9 = 4.5e3 varaux10 = 5e3 varaux11 = 5.5e3 varaux12 = 6e3 varaux13 = 6.5e3 varaux14 = 7e3 varaux15 = 7.5e3 varaux16 = 8e3 varaux17 = 8.5e3 varaux18 = 9e3 varaux19 = 9.5e3 varaux20 = 1.e4 varaux21 = 1.1e4 varaux22 = 1.2e4 varaux23 = 1.3e4 varaux24 = 1.4e4 varaux25 = 1.5e4 varaux26 = 1.7e4 varaux27 = 2.e4 varaux28 = 2.2e4 IF(tdimen/3.154e7.GE.0.AND.salida.EQ.0) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX1.AND.salida.EQ.1) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX2.AND.salida.eq.2) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX3.AND.salida.eq.3) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX4.AND.salida.eq.4) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX5.AND.salida.eq.5) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX6.AND.salida.eq.6) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX7.AND.salida.eq.7) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX8.AND.salida.eq.8) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX9.AND.salida.eq.9) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX10.AND.salida.eq.10) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX11.AND.salida.eq.11) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX12.AND.salida.eq.12) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX13.AND.salida.eq.13) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX14.AND.salida.eq.14) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX15.AND.salida.eq.15) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX16.AND.salida.eq.16) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX17.AND.salida.eq.17) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX18.AND.salida.eq.18) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX19.AND.salida.eq.19) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX20.AND.salida.eq.20) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX21.AND.salida.eq.21) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX22.AND.salida.eq.22) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX23.AND.salida.eq.23) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX24.AND.salida.eq.24) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX25.AND.salida.eq.25) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX26.AND.salida.eq.26) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX27.AND.salida.eq.27) imprimir1 = 1 IF(tdimen/3.154e7.GE.VARAUX28.AND.salida.eq.28) imprimir1 = 1 if(imprimir1.eq.1)then imprimir1 = 0 salida = salida + 1 imprime = imprime + 1 WRITE(1505,*) 'TITLE = "TodoAdim"' WRITE(1505,*) 'VARIABLES="X","X","U","V","T","C"' WRITE(1505,*) 'TEXT T="dimensionless time=',Tiempo & ,'" ZN=',IMPRIME WRITE(1505,*) 'ZONE T="',TIEMPO,'" I=',L1,', J=',M1 DO 850 J=1,M1 DO 850 I=1,L1 t(1,j) = t(2,j) t(l1,j) = t(l2,j) t(i,1) = t(i,2) c(i,1) = c(i,2) WRITE(1505,*) X(I),Y(J),U(I,J),V(I,J),T(I,J),C(I,J) 850 CONTINUE WRITE(1506,*) 'TITLE = "TodoDim"' WRITE(1506,*) 'VARIABLES="x (m)","y (m)","u (m/s)","v (m/s)" &,"T (K)","c (w/w)"' WRITE(1506,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1506,*) 'ZONE T="',(tdimen/3.154e7)/1000,'" I=',L1,', J=',M1 DO 840 J=1,M1 DO 840 I=1,L1 t(1,j) = t(2,j) t(l1,j) = t(l2,j) t(i,1) = t(i,2) c(i,1) = c(i,2) WRITE(1506,*) X(I)*LongC,Y(J)*LongC-6000,U0*U(I,J),U0*V(I,J), 1(((TINIC-TO)*T(I,J))+TO), 2(((CINIC-CO)*C(I,J))+CO) 840 CONTINUE WRITE(70,*) 'TITLE = "ConcDim"' WRITE(70,*) 'VARIABLES="x (m)","y (m)","c (w/w)"' WRITE(70,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(70,*) 'ZONE T="',(tdimen/3.154e7)/1000,'" I=',L1,', J=',M1 DO 700 J=1,M1 DO 700 I=1,L1 c(i,1) = c(i,2) WRITE(70,*) X(I)*LongC,Y(J)*LongC-6000,(((CINIC-CO)*C(I,J))+CO) 700 CONTINUE WRITE(71,*) 'TITLE = "TempDim"' WRITE(71,*) 'VARIABLES="x (m)","y (m)","T (K)"' WRITE(71,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(71,*) 'ZONE T="',(tdimen/3.154e7)/1000,'" I=',L1,', J=',M1 DO 701 J=1,M1 DO 701 I=1,L1 t(1,j) = t(2,j) t(l1,j) = t(l2,j) t(i,1) = t(i,2) WRITE(71,*) X(I)*LongC,Y(J)*LongC-6000,(((TINIC-TO)*T(I,J))+TO) 701 CONTINUE WRITE(72,*) 'TITLE = "VelocitiesDim"' WRITE(72,*) 'VARIABLES="x (m)","y (m)","u (m/s)","v (m/s)"' WRITE(72,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(72,*) 'ZONE T="',(tdimen/3.154e7)/1000,'" I=',L1,', J=',M1 DO 702 J=1,M1 DO 702 I=1,L1 WRITE(72,*) X(I)*LongC,y(j)*LongC-6000,U0*U(I,J),U0*V(I,J) 702 CONTINUE WRITE(1507,*) 'TITLE = "Velocity U Profile x"' WRITE(1507,*) 'VARIABLES="x (m)","u (m/s)"' WRITE(1507,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1507,*) 'ZONE t="',(tdimen/3.154e7)/1000,'" I=',L1 DO 703 I=1,L1 WRITE(1507,*) X(I)*LongC,U0*U(I,yprof1) 703 CONTINUE WRITE(1607,*) 'TITLE = "Velocity v Profile x"' WRITE(1607,*) 'VARIABLES="x (m)","v (m/s)"' WRITE(1607,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1607,*) 'ZONE t="',(tdimen/3.154e7)/1000,'" I=',L1 DO I=1,L1 WRITE(1607,*) X(I)*LongC,U0*V(i,yprof1) enddo WRITE(1508,*) 'TITLE = "Temp Profile x"' WRITE(1508,*) 'VARIABLES="x (m)","T (K)"' WRITE(1508,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1508,*) 'ZONE T="',(tdimen/3.154e7)/1000,'" I=',L1 DO 704 i=1,l1 WRITE(1508,*) X(I)*LongC,((TINIC-TO)*T(I,yprof1))+TO 704 CONTINUE WRITE(1608,*) 'TITLE = "Conc Profile x"' WRITE(1608,*) 'VARIABLES="x (m)","c (w/w)"' WRITE(1608,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1608,*) 'ZONE T="',(tdimen/3.154e7)/1000,'" I=',L1 DO i=1,l1 WRITE(1608,*) X(I)*LongC,((CINIC-CO)*C(I,xprof1))+CO enddo WRITE(1509,*) 'TITLE = "Velocities U Profile y"' WRITE(1509,*) 'VARIABLES="y (m)","u (m/s)"' WRITE(1509,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1509,*) 'ZONE t="',(tdimen/3.154e7)/1000,'" j=',m1 DO 705 j=1,m1 WRITE(1509,*) y(j)*LongC-6000,U0*U(xprof1,j) 705 CONTINUE WRITE(1609,*) 'TITLE = "Velocities V Profile y"' WRITE(1609,*) 'VARIABLES="y (m)","v (m/s)"' WRITE(1609,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1609,*) 'ZONE t="',(tdimen/3.154e7)/1000,'" j=',m1 DO j=1,m1 WRITE(1609,*) y(j)*LongC-6000,U0*V(xprof1,j) enddo WRITE(1510,*) 'TITLE = "Temp Profile y"' WRITE(1510,*) 'VARIABLES="y (m)","T (K)"' WRITE(1510,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1510,*) 'ZONE T="',(tdimen/3.154e7)/1000,'" j=',m1 DO 706 J=1,M1 WRITE(1510,*) y(j)*LongC-6000,((TINIC-TO)*T(xprof1,j))+TO 706 CONTINUE WRITE(1610,*) 'TITLE = "Conc Profile y"' WRITE(1610,*) 'VARIABLES="y (m)","c (w/w)"' WRITE(1610,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1610,*) 'ZONE T="',(tdimen/3.154e7)/1000,'" j=',m1 DO J=1,M1 WRITE(1610,*) y(j)*LongC-6000,((CINIC-CO)*C(xprof1,j))+CO enddo conti = 0 TPROM = 0 cprom = 0 uprom = 0 vprom = 0 vresprom =0 etaprom = 0. DO J=1,M1 DO I=1,L1 ! if(partes(i,j).eq.0.or.partes(i,j).eq.2)then if(camara(i,j).eq.0)then TPROM = T(I,J)+TPROM cprom = c(I,J)+cPROM uprom = u(I,J)+uPROM vprom = v(I,J)+vPROM etaprom = etai(I,J)+etaPROM Vresprom=Vresprom+SQRT(V(I,J)**2+U(I,J)**2) CONTI = CONTI + 1 !NUMERO DE VOLUMENES DE CONTROL CAMARA endif enddo enddo tprom = tprom/conti cprom = cprom/conti uprom = uprom/conti vprom = vprom/conti etaprom = etaprom/conti vresprom = vresprom/conti WRITE(73,*) (tdimen/3.154e7)/1000,((TINIC-TO)*tprom)+TO,tprom WRITE(74,*) (tdimen/3.154e7)/1000,((CINIC-CO)*cprom)+CO,cprom WRITE(75,*) (tdimen/3.154e7)/1000,U0*uprom,uprom WRITE(76,*) (tdimen/3.154e7)/1000,U0*vprom,vprom WRITE(77,*) (tdimen/3.154e7)/1000,etaprom WRITE(88,*) (tdimen/3.154e7)/1000,U0*vresprom,vresprom ENDIF ! write(*,*)tprom,MaxTemp,(((TINIC-TO)*MaxTemp)+TO) ! pause ! IF(TIEMPO.GE.TFIN+DTPO.or.tprom.le.0.958)THEN IF(TIEMPO.GE.TFIN+DTPO &.or.((TINIC-TO)*T(xprof1,yprof1))+TO.le.1448)THEN conti = 0 TPROM = 0 cprom = 0 uprom = 0 vprom = 0 vresprom =0 etaprom = 0. DO J=1,M1 DO I=1,L1 ! if(partes(i,j).eq.0.or.partes(i,j).eq.2)then if(camara(i,j).eq.0)then TPROM = T(I,J)+TPROM cprom = c(I,J)+cPROM uprom = u(I,J)+uPROM vprom = v(I,J)+vPROM etaprom = etai(I,J)+etaPROM Vresprom=Vresprom+SQRT(V(I,J)**2+U(I,J)**2) CONTI = CONTI + 1 !NUMERO DE VOLUMENES DE CONTROL CAMARA endif enddo enddo tprom = tprom/conti cprom = cprom/conti uprom = uprom/conti vprom = vprom/conti etaprom = etaprom/conti vresprom = vresprom/conti WRITE(73,*) (tdimen/3.154e7)/1000,((TINIC-TO)*tprom)+TO,tprom WRITE(74,*) (tdimen/3.154e7)/1000,((CINIC-CO)*cprom)+CO,cprom WRITE(75,*) (tdimen/3.154e7)/1000,U0*uprom,uprom WRITE(76,*) (tdimen/3.154e7)/1000,U0*vprom,vprom WRITE(77,*) (tdimen/3.154e7)/1000,etaprom WRITE(88,*) (tdimen/3.154e7)/1000,U0*vresprom,vresprom WRITE(1504,*) 'TITLE = "tiempo final dimensional"' WRITE(1504,*) 'VARIABLES="x (m)","y (m)","u (m/s)","v (m/s)" & ,"T (k)","c (w/w)"' WRITE(1504,*) 'TEXT T="ky=',(tdimen/3.154e7)/1000 & ,'" ZN=',IMPRIME WRITE(1504,*) 'ZONE T="',(tdimen/3.154e7)/1000,' & " I=',L1,', J=',M1 DO J=1,M1 DO I=1,L1 t(1,j) = t(2,j) t(l1,j) = t(l2,j) t(i,1) = t(i,2) c(i,1) = c(i,2) WRITE(1504,*) X(I)*LongC,Y(J)*LongC-6000,U0*U(I,J),U0*V(I,J), 1 (((TINIC-TO)*T(I,J))+TO), 2 (((CINIC-CO)*C(I,J))+CO) enddo enddo DTIEMPO=ETIME(TA)*(tiempo * LongC)/U0 WRITE(*,*) 'El programa funciono por:',DTIEMPO,'CPU segundos.' WRITE(*,*) 'DONDE:' WRITE(*,*) 'Tpo de usuario',TA(1),'segundos' WRITE(*,*) 'Tpo de sistema:',TA(2),'segundos' WRITE(*,*) 'Tpo usuario',TA(1),'anos' WRITE(*,*) 'Tpo sistema',TA(2),'anos' WRITE(*,*) 'La cantidad total de iteraciones fue:',ITER WRITE(*,*) 'tpo adim:',tiempo,'tpo dim:',tdimen, & 'en años:',tdimen/3.154e7 WRITE(101,*) 'El programa funciono por:',DTIEMPO,'CPU segundos.' WRITE(101,*) 'DONDE:' WRITE(101,*) 'Tpo usuario',TA(1),'segundos' WRITE(101,*) 'Tpo sistema:',TA(2),'segundos' WRITE(101,*) 'Tpo usrio',TA(1),'años' WRITE(101,*) 'Tpo sist',TA(2),'años' WRITE(101,*) 'La cantidad total de iteraciones fue:',ITER WRITE(101,*) 'tpo adim:',tiempo,'tpo dim:',tdimen, & 'en años:',tdimen/3.154e7 WRITE(101,*) 'Volumen Camara en m2=',volumen & ,'en Km2=',Volumen/1e6, & 'LCarac m=',LongC PAUSE STOP ENDIF DO 470 I=1,L1 DO 470 J=1,M1 UANT(I,J)=U(I,J) VANT(I,J)=V(I,J) TANT(I,J)=T(I,J) cANT(I,J)=c(I,J) 470 CONTINUE !******SE CALCULA LA TEMPERATURA PROMEDIO DE LA CÁMARA ITERN=0. TIEMPO=TIEMPO+DTPO ENDIF ITERN=ITERN+1 RETURN C**************************************************************** C**************************************************************** ENTRY GAMSOR C**********calculo liquido interior no newtoniano ********************* ! etai = 1e30 IF(ITERN.GT.5.AND.AN.NE.1.)THEN DO 502 I=2,L2 DO 502 J=2,M2 IF(partes(I,J).EQ.0.or.partes(i,j).eq.2)THEN *******CALCULO DE dU/dY************************************************ UP=0.5/YDIF(J+1)*(U(I,J+1)*YCV(J+1)+U(I,J)*YCV(J)) UM=0.5/YDIF(J)*(U(I,J)*YCV(J)+U(I,J-1)*YCV(J-1)) IF(J.EQ.2) UM=U(I,J-1) IF(J.EQ.M2) UP=U(I,J+1) DUDYS=(UP-UM)/(YCV(J)+SMALL) UP=0.5/YDIF(J+1)*(U(I+1,J+1)*YCV(J+1)+U(I+1,J)*YCV(J)) UM=0.5/YDIF(J)*(U(I+1,J)*YCV(J)+U(I+1,J-1)*YCV(J-1)) IF(J.EQ.2) UM=U(I,J-1) IF(J.EQ.M2) UP=U(I,J+1) DUDYN=(UP-UM)/(YCV(J)+SMALL) DUDY=0.5*(DUDYS+DUDYN) *******CALCULO DE dV/dX************************************************ VP=0.5/XDIF(I+1)*(V(I+1,J)*XCV(I+1)+V(I,J)*XCV(I)) VM=0.5/XDIF(I)*(V(I,J)*XCV(I)+V(I-1,J)*XCV(I-1)) IF(I.EQ.2) VM=V(I-1,J) IF(I.EQ.L2) VP=V(I+1,J) DVDXW=(VP-VM)/(XCV(I)+SMALL) VP=0.5/XDIF(I+1)*(V(I+1,J+1)*XCV(I+1)+V(I,J+1)*XCV(I)) VM=0.5/XDIF(I)*(V(I,J+1)*XCV(I)+V(I-1,J+1)*XCV(I-1)) IF(I.EQ.2) VM=V(I-1,J) IF(I.EQ.L2) VP=V(I+1,J) DVDXE=(VP-VM)/(XCV(I)+SMALL) DVDX=0.5*(DVDXW+DVDXE) c*******CALCULO DE dV/dY************************************************ DVDY=(V(I,J+1)-V(I,J))/YCV(J) IF(J.EQ.M2)DVDY=(V(I,J+1)-V(I,J))/(0.5*YCV(J)) c*******CALCULO DE dU/dX************************************************ DUDX=(U(I+1,J)-U(I,J))/XCV(I) IF(I.EQ.L2)DUDX=(U(I+1,J)-U(I,J))/(0.5*XCV(I)) c*********************************************************************** GAMAP=2*(DVDY**2+DUDX**2)+(DVDX+DUDY)**2 F(I,J,6)=SQRT(GAMAP)*U0/LongC !shear rate o tasa de deformacion, y punto adim IF(F(I,J,6).LT.0.01) F(I,J,6)=0.01 !minima tasa de deformación dim c****aplica los esfuerzos de corte a la viscosidad aparente power law model**** if(abs(((TINIC-TO)*T(I,J))+TO).lt.1448)then ! etai(i,j) = 0. ETAI(I,J)=1e+30 !/ETAREF partes(i,j) = 1 goto 321 endif if(abs(((TINIC-TO)*T(I,J))+TO).gt.1498)then ! etai(i,j) = 0. ETAI(I,J)=0.1/ETAREF ! partes(i,j) = 1 goto 321 endif dE(I,J)=YY0+AA0*abs(((TINIC-TO)*T(I,J))+TO) !energía activación dim F(I,J,7)=EXP((AN*dE(I,J))/(CTEGAS*abs((T(I,J) !exponente de formula para densidad aparente & *(TINIC-TO)+TO)))) !densidad aparente dim ETAI(I,J)=(ETA0*F(I,J,6)**(AN-1)*F(I,J,7))/ETAREF !dens aparente adim. eta0 es el indice de consistencia 321 continue endif 502 continue endif !original ! FuenteVel= sqrt(ra0/pr0) !*sqrt(LongCL/LongC) !U0**2/LongC ! FuenteVelFlota= sqrt(Ra0/pr0) !sqrt(Ra0/pr0)*(LongC/LongCL) ! GamModelV = 1. !sqrt(longC/LongCL)/sqrt(ra0/pr0) ! GamModelTS = 1. !sqrt(LongC/LongCL)/(sqrt(ra0*pr0)*alfastar) ! GamModelT = 1. !sqrt(LongC/LongCL)/sqrt(ra0*pr0) ! FuenteT = sqrt(ra0*pr0) !sqrt(Ra0*Pr0)*sqrt(LongCL/LongC) ! FuenteTS = sqrt(Ra0*Pr0)*alfastar !sqrt(Ra0*Pr0)*alfastar*sqrt(LongCL/LongC) ! GamModelC = 1. !sqrt(longC/LongCL)/sqrt(ra0/pr0 ! FuenteC = sqrt(Ra0*Pr0)*Le !sqrt(Ra0*Pr0)*Le*sqrt(LongCL/LongC) ! FuenteCS = sqrt(Ra0*Pr0)*Le !sqrt(Ra0*Pr0)*Le*sqrt(LongCL/LongC) !nuevo modificado 1 (distribuido en el lado del tiempo y en gam GamModelV = sqrt(longC/LongCL) FuenteVel= sqrt(ra0/pr0) FuenteVelFlota= sqrt(Ra0/pr0)*(LongC/LongCL)*(1./deltaTad) GamModelTS = sqrt(longC/LongCL) GamModelT = sqrt(longC/LongCL) FuenteT = sqrt(ra0*pr0) FuenteTS = sqrt(Ra0*Pr0)*alfastar GamModelC = sqrt(longC/LongCL) FuenteC = sqrt(Ra0*Pr0)*Le FuenteCS = sqrt(Ra0*Pr0)*Le !nuevo modificado 2 (todo al lado del gama) ! GamModelV = sqrt(longC/LongCL)/sqrt(ra0/pr0) ! FuenteVel= 1. ! FuenteVelFlota= (LongC/LongCL) ! GamModelTS = sqrt(LongC/LongCL)/(sqrt(ra0*pr0)*alfastar) ! GamModelT = sqrt(LongC/LongCL)/sqrt(ra0*pr0) ! FuenteT = 1. ! FuenteTS = 1. ! GamModelC = sqrt(longC/LongCL)/sqrt(ra0/pr0) ! FuenteC = 1. ! FuenteCS = 1. DO 530 I=2,L2 DO 530 J=2,M2 C COMPONENTE U ********************** IF(NF.EQ.1) THEN CON(I,J) = UANT(I,J)/DTPO AP(I,J) = -1./DTPO if(partes(i,j).eq.1)gam(i,j) = 1E30 if(partes(i,j).eq.0.or.partes(i,j).eq.2)then gam(i,1) = etai(i,1)*GamModelV gam(i,j) = etai(i,j)*GamModelV CON(I,J) = UANT(I,J)*FuenteVel/DTPO AP(I,J) = -1.*FuenteVel/DTPO endif ! gam(i,1)=1e30 ENDIF C COMPONENTE V ********************** IF(NF.EQ.2) THEN CON(I,J) = (VANT(I,J)/DTpo) AP(I,J) = -1./DTPO if(partes(i,j).eq.1)gam(i,j) = 1E30 if(partes(i,j).eq.0.or.partes(i,j).eq.2)then gam(i,1) = etai(i,1)*GamModelV gam(i,j) = etai(i,j)*GamModelV IF (i.EQ.1) TM = T(I,J) IF (i.NE.1) TM = 0.5*(T(I,J)+T(I,J-1)) IF (j.GT.2) CON(I,J) = VANT(I,J)*FuenteVel/DTpo & - FuenteVelFlota*TM AP(I,J) = -1.*FuenteVel/DTPO endif ! gam(i,1)=1e30 ENDIF 530 CONTINUE C TEMPERATURA ********************** DO 540 I=1,L1 DO 540 J=1,M1 IF(NF.EQ.4)THEN CON(I,J) = TANT(I,J)*FuenteTS/DTPO AP(I,J) = - FuenteTS/DTPO if(partes(i,j).eq.1)gam(i,j) = 1.*GamModelTS !solido if(partes(i,j).eq.0.or.partes(i,j).eq.2)then !liquido gam(i,j) = 1.* GamModelT CON(I,J) = TANT(I,J)*FuenteT/DTPO AP(I,J) = - FuenteT/DTPO endif gam(1,j) = 0. gam(l1,j)= 0. ! gam(i,1)= 0. ENDIF 540 CONTINUE C CONCENTRACION ********************** DO 550 I=1,L1 DO 550 J=1,M1 IF(NF.EQ.5)THEN CON(I,J) = cANT(I,J)*FuenteC/DTPO AP(I,J) = - 1.*FuenteC/DTPO if(partes(i,j).eq.1)gam(i,j) = 0 !1.*GamModelC if(partes(i,j).eq.0.or.partes(i,j).eq.2)then gam(i,j) = 1.*GamModelC CON(I,J) = cANT(I,J)*FuenteC/DTPO AP(I,J) = - FuenteC/DTPO endif gam(1,j) = 0. gam(l1,j)= 0. ! gam(i,1)= 0. ENDIF 550 CONTINUE RETURN END !********************************* !SUBPROGRAMA GENERADOR DE MALLA ! VARIABLE A TRAMOS !******************************** SUBROUTINE MALLA_VAT !******************************** PARAMETER (NOD=500) !NUMERO MAXIMO DE NODOS COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, 1IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, 2IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCO COMMON F(NOD,NOD,10),P(NOD,NOD),RHO(NOD,NOD),GAM(NOD,NOD), 1CON(NOD,NOD), 1AIP(NOD,NOD),AIM(NOD,NOD),AJP(NOD,NOD),AJM(NOD,NOD),AP(NOD,NOD), 2X(NOD),XU(NOD),XDIF(NOD),XCV(NOD),XCVS(NOD), 3Y(NOD),YV(NOD),YDIF(NOD),YCV(NOD),YCVS(NOD), 4YCVR(NOD),YCVRS(NOD),ARX(NOD),ARXJ(NOD),ARXJP(NOD), 5R(NOD),RMN(NOD),SX(NOD),SXMN(NOD),XCVI(NOD),XCVIP(NOD), 7YLS(10),M1S(10),M1ST(10),XLS(10),L1S(10),L1ST(10),NSX,NSY,RAX(10) 8,RAY(10) COMMON/MVAT/NTYV,NTXV,VTYV(30),VTXV(30),YLTV(30),XLTV(30),ITRAM, 1COORD,RX(30),RY(30) !*************************************************************** XL=0. DO I=1,INT(NTXV) XL=XL+XLTV(I) ENDDO L1=2 DO I=1,INT(NTXV) L1=L1+INT(VTXV(I)) ENDDO !*************************************************************** K=1 XNN=VTXV(K) XNN_1=0 DO I=0,L1-3 IF(RX(K).NE.1)THEN A=XLTV(K)*(RX(K)-1.)/(RX(K)**(VTXV(K))-1.) XU(I+3)=A*RX(K)**(FLOAT(I)-XNN_1)+XU(I+2) ELSE DX=XLTV(K)/VTXV(K) XU(I+3)=XU(I+2)+DX ENDIF IF(FLOAT(I).EQ.(XNN-1.))THEN XNN_1=XNN_1+VTXV(K) K=K+1 XNN=XNN+VTXV(K) ENDIF ENDDO !*************************************************************** YL=0. DO I=1,INT(NTYV) YL=YL+YLTV(I) ENDDO M1=2 DO I=1,INT(NTYV) M1=M1+INT(VTYV(I)) ENDDO !*************************************************************** K=1 XNN=VTYV(K) XNN_1=0 DO I=0,M1-3 IF(RY(K).NE.1)THEN B=YLTV(K)*(RY(K)-1.)/(RY(K)**(VTYV(K))-1.) YV(I+3)=B*RY(K)**(FLOAT(I)-XNN_1)+YV(I+2) ELSE DY=YLTV(K)/VTYV(K) YV(I+3)=YV(I+2)+DY ENDIF IF(FLOAT(I).EQ.(XNN-1.))THEN XNN_1=XNN_1+VTYV(K) K=K+1 XNN=XNN+VTYV(K) ENDIF ENDDO RETURN END C LINE0010 C***********************************************************************LINE0020 C LINE0030 C LINE0040 C-------------------------- METODO SIMPLE-----------------------------LINE0050 C LINE0060 C_______________________ VOLUMENES DE CONTROL---------------------------LINE0070 C LINE0080 C***********************************************************************LINE0090 C LINE0100 c INTEGER H,HI,HF,M,MI,MF,S,SI,SF,C,CI,CF LOGICAL LSTOP LINE0110 COMMON/CNTL/LSTOP LINE0120 C character*12 tinicio,tfinal c CALL TIME (HI,MI,SI,CI) C call time(tinicio) CALL OPENF CALL GRID LINE0130 CALL SETUP1 LINE0140 CALL START LINE0150 10 CALL DENSE LINE0160 CALL BOUND LINE0170 CALL OUTPUT LINE0180 IF(LSTOP) THEN STOP END IF LINE0190 CALL SETUP2 LINE0200 GO TO 10 END SUBROUTINE DIFLOW LINE0230 COMMON/COEF/FLOW,DIFF,ACOF LINE0240 ACOF=DIFF LINE0250 IF(FLOW.EQ.0.) RETURN LINE0260 TEMP=DIFF-ABS(FLOW)*0.1 LINE0270 ACOF=0. LINE0280 IF(TEMP.LE.0.)RETURN LINE0290 TEMP=TEMP/DIFF LINE0300 ACOF=DIFF*TEMP**5 LINE0310 RETURN LINE0320 END LINE0330 SUBROUTINE SOLVE LINE0340 CHARACTER*10 TITLE LINE0350 LOGICAL LSOLVE,LPRINT,LBLK !,LSTOP LINE0360 PARAMETER (NOD=500) !NUMERO MAXIMO DE NODOS COMMON F(NOD,NOD,10),P(NOD,NOD),RHO(NOD,NOD),GAM(NOD,NOD), LINE0365 1 CON(NOD,NOD), LINE0370 2 AIP(NOD,NOD),AIM(NOD,NOD),AJP(NOD,NOD),AJM(NOD,NOD),AP(NOD,NOD), LINE0380 3 X(NOD),XU(NOD),XDIF(NOD),XCV(NOD),XCVS(NOD), LINE0390 4 Y(NOD),YV(NOD),YDIF(NOD),YCV(NOD),YCVS(NOD), LINE0400 5 YCVR(NOD),YCVRS(NOD),ARX(NOD),ARXJ(NOD),ARXJP(NOD), LINE0410 6 R(NOD),RMN(NOD),SX(NOD),SXMN(NOD),XCVI(NOD),XCVIP(NOD), LINE0420 7 YLS(10),M1S(10),M1ST(10),XLS(10),L1S(10),L1ST(10),NSX,NSY COMMON DU(NOD,NOD),DV(NOD,NOD),FV(NOD),FVP(NOD), LINE0430 1 FX(NOD),FXM(NOD),FY(NOD),FYM(NOD),PT(NOD),QT(NOD) LINE0440 COMMON/CHAR/TITLE(13) LINE0450 COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, LINE0460 1IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, LINE0470 2IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON LINE0480 ISTF=IST-1 LINE0490 JSTF=JST-1 LINE0500 IT1=L2+IST LINE0510 IT2=L3+IST LINE0520 JT1=M2+JST LINE0530 JT2=M3+JST LINE0540 DO 999 NT=1,NTIMES(NF) LINE0550 DO 999 N=NF,NF LINE0560 C---------------------------------------------------------------------- LINE0570 IF(.NOT.LBLK(NF)) GO TO 10 LINE0580 PT(ISTF)=0. LINE0590 QT(ISTF)=0. LINE0600 DO 11 I=IST,L2 LINE0610 BL=0. LINE0620 BLP=0. LINE0630 BLM=0. LINE0640 BLC=0. LINE0650 DO 12 J=JST,M2 LINE0660 BL=BL+AP(I,J) LINE0670 IF(J.NE.M2) BL=BL-AJP(I,J) LINE0680 IF(J.NE.JST) BL=BL-AJM(I,J) LINE0690 BLP=BLP+AIP(I,J) LINE0700 BLM=BLM+AIM(I,J) LINE0710 BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) LINE0720 1 +AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N) LINE0730 12 CONTINUE LINE0740 DENOM=BL-PT(I-1)*BLM LINE0750 IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E20 LINE0760 PT(I)=BLP/DENOM LINE0770 QT(I)=(BLC+BLM*QT(I-1))/DENOM LINE0780 11 CONTINUE LINE0790 BL=0. LINE0800 DO 13 II=IST,L2 LINE0810 I=IT1-II LINE0820 BL=BL*PT(I)+QT(I) LINE0830 DO 13 J=JST,M2 LINE0840 13 F(I,J,N)=F(I,J,N)+BL LINE0850 PT(JSTF)=0. LINE0860 QT(JSTF)=0. LINE0870 DO 21 J=JST,M2 LINE0880 BL=0. LINE0890 BLP=0. LINE0900 BLM=0. LINE0910 BLC=0. LINE0920 DO 22 I=IST,L2 LINE0930 BL=BL+AP(I,J) LINE0940 IF(I.NE.L2) BL=BL-AIP(I,J) LINE0950 IF(I.NE.IST) BL=BL-AIM(I,J) LINE0960 BLP=BLP+AJP(I,J) LINE0970 BLM=BLM+AJM(I,J) LINE0980 BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) LINE0990 1 +AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N) LINE1000 22 CONTINUE LINE1010 DENOM=BL-PT(J-1)*BLM LINE1020 IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E20 LINE1030 PT(J)=BLP/DENOM LINE1040 QT(J)=(BLC+BLM*QT(J-1))/DENOM LINE1050 21 CONTINUE LINE1060 BL=0. LINE1070 DO 23 JJ=JST,M2 LINE1080 J=JT1-JJ LINE1090 BL=BL*PT(J)+QT(J) LINE1100 DO 23 I=IST,L2 LINE1110 23 F(I,J,N)=F(I,J,N)+BL LINE1120 10 CONTINUE LINE1130 DO 90 J=JST,M2 LINE1140 PT(ISTF)=0. LINE1150 QT(ISTF)=F(ISTF,J,N) LINE1160 DO 70 I=IST,L2 LINE1170 DENOM=AP(I,J)-PT(I-1)*AIM(I,J) LINE1180 PT(I)=AIP(I,J)/DENOM LINE1190 TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N) LINE1200 QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM LINE1210 70 CONTINUE LINE1220 DO 80 II=IST,L2 LINE1230 I=IT1-II LINE1240 80 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I) LINE1250 90 CONTINUE LINE1260 DO 190 JJ=JST,M3 LINE1270 J=JT2-JJ LINE1280 PT(ISTF)=0. LINE1290 QT(ISTF)=F(ISTF,J,N) LINE1300 DO 170 I=IST,L2 LINE1310 DENOM=AP(I,J)-PT(I-1)*AIM(I,J) LINE1320 PT(I)=AIP(I,J)/DENOM LINE1330 TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N) LINE1340 QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM LINE1350 170 CONTINUE LINE1360 DO 180 II=IST,L2 LINE1370 I=IT1-II LINE1380 180 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I) LINE1390 190 CONTINUE LINE1400 DO 290 I=IST,L2 LINE1410 PT(JSTF)=0. LINE1420 QT(JSTF)=F(I,JSTF,N) LINE1430 DO 270 J=JST,M2 LINE1440 DENOM=AP(I,J)-PT(J-1)*AJM(I,J) LINE1450 PT(J)=AJP(I,J)/DENOM LINE1460 TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) LINE1470 QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM LINE1480 270 CONTINUE LINE1490 DO 280 JJ=JST,M2 LINE1500 J=JT1-JJ LINE1510 280 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J) LINE1520 290 CONTINUE LINE1530 DO 390 II=IST,L3 LINE1540 I=IT2-II LINE1550 PT(JSTF)=0. LINE1560 QT(JSTF)=F(I,JSTF,N) LINE1570 DO 370 J=JST,M2 LINE1580 DENOM=AP(I,J)-PT(J-1)*AJM(I,J) LINE1590 PT(J)=AJP(I,J)/DENOM LINE1600 TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) LINE1610 QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM LINE1620 370 CONTINUE LINE1630 DO 380 JJ=JST,M2 LINE1640 J=JT1-JJ LINE1650 380 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J) LINE1660 390 CONTINUE LINE1670 999 CONTINUE LINE1680 DO 400 J=2,M2 LINE1690 DO 400 I=2,L2 LINE1700 CON(I,J)=0. LINE1710 AP(I,J)=0. LINE1720 400 CONTINUE LINE1730 RETURN LINE1740 END LINE1750 SUBROUTINE SETUP LINE1760 CHARACTER*10 TITLE LINE1770 LOGICAL LSOLVE,LPRINT,LBLK,LSTOP LINE1780 PARAMETER (NOD=500) !NUMERO MAXIMO DE NODOS COMMON F(NOD,NOD,10),P(NOD,NOD),RHO(NOD,NOD),GAM(NOD,NOD), LINE1785 1 CON(NOD,NOD), LINE1790 2 AIP(NOD,NOD),AIM(NOD,NOD),AJP(NOD,NOD),AJM(NOD,NOD),AP(NOD,NOD), LINE1800 3 X(NOD),XU(NOD),XDIF(NOD),XCV(NOD),XCVS(NOD), LINE1810 4 Y(NOD),YV(NOD),YDIF(NOD),YCV(NOD),YCVS(NOD), LINE1820 5 YCVR(NOD),YCVRS(NOD),ARX(NOD),ARXJ(NOD),ARXJP(NOD), LINE1830 6 R(NOD),RMN(NOD),SX(NOD),SXMN(NOD),XCVI(NOD),XCVIP(NOD), LINE1840 7 YLS(10),M1S(10),M1ST(10),XLS(10),L1S(10),L1ST(10),NSX,NSY COMMON DU(NOD,NOD),DV(NOD,NOD),FV(NOD),FVP(NOD), LINE1850 1 FX(NOD),FXM(NOD),FY(NOD),FYM(NOD),PT(NOD),QT(NOD) LINE1860 COMMON/CHAR/TITLE(13) LINE1870 COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, LINE1880 1IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, LINE1890 2IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON LINE1900 COMMON/CNTL/LSTOP LINE1910 COMMON/SORC/SMAX,SSUM LINE1920 COMMON/COEF/FLOW,DIFF,ACOF LINE1930 REAL U(NOD,NOD),V(NOD,NOD),PC(NOD,NOD) LINE1940 C DIMENSION U(NOD,NOD),V(NOD,NOD),PC(NOD,NOD) LINE1940 EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1)) LINE1950 101 FORMAT(//,3X,'COMPUTATION IN CARTESIAN COORDINATES') LINE1960 2 FORMAT(//,3X,'COMPUTATION FOR AXISYMMETRIC SITUATION') LINE1970 3 FORMAT(//,3X,'COMPUTATION IN POLAR COORDINATES') LINE1980 4 FORMAT(2X,40(1H*),//) SI01990 ENTRY SETUP1 LINE2000 L2=L1-1 LINE2010 L3=L2-1 LINE2020 M2=M1-1 LINE2030 M3=M2-1 LINE2040 X(1)=XU(2) LINE2050 DO 5 I=2,L2 LINE2060 5 X(I)=0.5*(XU(I+1)+XU(I)) LINE2070 X(L1)=XU(L1) LINE2080 Y(1)=YV(2) LINE2090 DO 10 J=2,M2 LINE2100 10 Y(J)=0.5*(YV(J+1)+YV(J)) LINE2110 Y(M1)=YV(M1) LINE2120 DO 15 I=2,L1 LINE2130 15 XDIF(I)=X(I)-X(I-1) LINE2140 DO 18 I=2,L2 LINE2150 18 XCV(I)=XU(I+1)-XU(I) LINE2160 DO 20 I=3,L2 LINE2170 20 XCVS(I)=XDIF(I) LINE2180 XCVS(3)=XCVS(3)+XDIF(2) LINE2190 XCVS(L2)=XCVS(L2)+XDIF(L1) LINE2200 DO 22 I=3,L3 LINE2210 XCVI(I)=0.5*XCV(I) LINE2220 22 XCVIP(I)=XCVI(I) LINE2230 XCVIP(2)=XCV(2) LINE2240 XCVI(L2)=XCV(L2) LINE2250 DO 35 J=2,M1 LINE2260 35 YDIF(J)=Y(J)-Y(J-1) LINE2270 DO 40 J=2,M2 LINE2280 40 YCV(J)=YV(J+1)-YV(J) LINE2290 DO 45 J=3,M2 LINE2300 45 YCVS(J)=YDIF(J) LINE2310 YCVS(3)=YCVS(3)+YDIF(2) LINE2320 YCVS(M2)=YCVS(M2)+YDIF(M1) LINE2330 IF(MODE.NE.1) GO TO 55 LINE2340 DO 52 J=1,M1 LINE2350 RMN(J)=1.0 LINE2360 52 R(J)=1.0 LINE2370 GO TO 56 LINE2380 55 DO 50 J=2,M1 LINE2390 50 R(J)=R(J-1)+YDIF(J) LINE2400 RMN(2)=R(1) LINE2410 DO 60 J=3,M2 LINE2420 60 RMN(J)=RMN(J-1)+YCV(J-1) LINE2430 RMN(M1)=R(M1) LINE2440 56 CONTINUE LINE2450 DO 57 J=1,M1 LINE2460 SX(J)=1. LINE2470 SXMN(J)=1. LINE2480 IF(MODE.NE.3) GO TO 57 LINE2490 SX(J)=R(J) LINE2500 IF(J.NE.1) SXMN(J)=RMN(J) LINE2510 57 CONTINUE LINE2520 DO 62 J=2,M2 LINE2530 YCVR(J)=R(J)*YCV(J) LINE2540 ARX(J)=YCVR(J) LINE2550 IF(MODE.NE.3) GO TO 62 LINE2560 ARX(J)=YCV(J) LINE2570 62 CONTINUE LINE2580 DO 64 J=4,M3 LINE2590 64 YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J) LINE2600 YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3) LINE2610 YCVRS(M2)=0.5*(R(M1)+R(M3))*YCVS(M2) LINE2620 IF(MODE.NE.2) GO TO 67 LINE2630 DO 65 J=3,M3 LINE2640 ARXJ(J)=0.25*(1.+RMN(J)/R(J))*ARX(J) LINE2650 65 ARXJP(J)=ARX(J)-ARXJ(J) LINE2660 GO TO 68 LINE2670 67 DO 66 J=3,M3 LINE2680 ARXJ(J)=0.5*ARX(J) LINE2690 66 ARXJP(J)=ARXJ(J) LINE2700 68 ARXJP(2)=ARX(2) LINE2710 ARXJ(M2)=ARX(M2) LINE2720 DO 70 J=3,M3 LINE2730 FV(J)=ARXJP(J)/ARX(J) LINE2740 70 FVP(J)=1.-FV(J) LINE2750 DO 85 I=3,L2 LINE2760 FX(I)=0.5*XCV(I-1)/XDIF(I) LINE2770 85 FXM(I)=1.-FX(I) LINE2780 FX(2)=0. LINE2790 FXM(2)=1. LINE2800 FX(L1)=1. LINE2810 FXM(L1)=0. LINE2820 DO 90 J=3,M2 LINE2830 FY(J)=0.5*YCV(J-1)/YDIF(J) LINE2840 90 FYM(J)=1.-FY(J) LINE2850 FY(2)=0. LINE2860 FYM(2)=1. LINE2870 FY(M1)=1. LINE2880 FYM(M1)=0. LINE2890 CCON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE LINE2900 DO 95 J=1,M1 LINE2910 DO 95 I=1,L1 LINE2920 PC(I,J)=0. LINE2930 U(I,J)=0. LINE2940 V(I,J)=0. LINE2950 CON(I,J)=0. LINE2960 AP(I,J)=0. LINE2970 RHO(I,J)=RHOCON LINE2980 P(I,J)=0. LINE2990 95 CONTINUE LINE3000 IF (MODE.EQ.1) then ! write(*,101) ! write(1,*)' Computation in cartesian coordinates' LINE3010 end if IF (MODE.EQ.2) then c PRINT 2 c write(1,*)'coord.cilindricas' LINE3020 end if IF (MODE.EQ.3) then c PRINT 3 c write(1,3) LINE3030 end if c PRINT 4 c write(1,4) LINE3040 RETURN LINE3050 ENTRY SETUP2 LINE3060 CCOEFFICIENTS FOR THE U EQUATION LINE3070 NF=1 LINE3080 IF(.NOT.LSOLVE(NF)) GO TO 100 LINE3090 IST=3 LINE3100 JST=2 LINE3110 CALL GAMSOR LINE3120 REL=1.-RELAX(NF) LINE3130 DO 102 I=3,L2 LINE3140 FL=XCVI(I)*V(I,2)*RHO(I,1) LINE3150 FLM=XCVIP(I-1)*V(I-1,2)*RHO(I-1,1) LINE3160 FLOW=R(1)*(FL+FLM) LINE3170 DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2) LINE3180 CALL DIFLOW LINE3190 102 AJM(I,2)=ACOF+AMAX1(0.,FLOW) LINE3200 DO 103 J=2,M2 LINE3210 FLOW=ARX(J)*U(2,J)*RHO(1,J) LINE3220 DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J)) LINE3230 CALL DIFLOW LINE3240 AIM(3,J)=ACOF+AMAX1(0.,FLOW) LINE3250 DO 103 I=3,L2 LINE3260 IF(I.EQ.L2) GO TO 104 LINE3270 FL=U(I,J)*(FX(I)*RHO(I,J)+FXM(I)*RHO(I-1,J)) LINE3280 FLP=U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) LINE3290 FLOW=ARX(J)*0.5*(FL+FLP) LINE3300 DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J)) LINE3310 GO TO 105 LINE3320 104 FLOW=ARX(J)*U(L1,J)*RHO(L1,J) LINE3330 DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J)) LINE3340 105 CALL DIFLOW LINE3350 AIM(I+1,J)=ACOF+AMAX1(0.,FLOW) LINE3360 AIP(I,J)=AIM(I+1,J)-FLOW LINE3370 IF(J.EQ.M2) GO TO 106 LINE3380 FL=XCVI(I)*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J)) LINE3390 FLM=XCVIP(I-1)*V(I-1,J+1)*(FY(J+1)*RHO(I-1,J+1)+FYM(J+1)* LINE3400 1 RHO(I-1,J)) LINE3410 GM=GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+ LINE3420 1 1.0E-30)*XCVI(I) LINE3430 GMM=GAM(I-1,J)*GAM(I-1,J+1)/(YCV(J)*GAM(I-1,J+1)+YCV(J+1)* LINE3440 1 GAM(I-1,J)+1.E-30)*XCVIP(I-1) LINE3450 DIFF=RMN(J+1)*2.*(GM+GMM) LINE3460 GO TO 107 LINE3470 106 FL=XCVI(I)*V(I,M1)*RHO(I,M1) LINE3480 FLM=XCVIP(I-1)*V(I-1,M1)*RHO(I-1,M1) LINE3490 DIFF=R(M1)*(XCVI(I)*GAM(I,M1)+XCVIP(I-1)*GAM(I-1,M1))/YDIF(M1) LINE3500 107 FLOW=RMN(J+1)*(FL+FLM) LINE3510 CALL DIFLOW LINE3520 AJM(I,J+1)=ACOF+AMAX1(0.,FLOW) LINE3530 AJP(I,J)=AJM(I,J+1)-FLOW LINE3540 VOL=YCVR(J)*XCVS(I) LINE3550 APT=(RHO(I,J)*XCVI(I)+RHO(I-1,J)*XCVIP(I-1)) LINE3560 1/(XCVS(I)*DT) LINE3570 AP(I,J)=AP(I,J)-APT LINE3580 CON(I,J)=CON(I,J)+APT*U(I,J) LINE3590 AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) LINE3600 1/RELAX(NF) LINE3610 CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*U(I,J) LINE3620 DU(I,J)=VOL/(XDIF(I)*SX(J)) LINE3630 CON(I,J)=CON(I,J)+DU(I,J)*(P(I-1,J)-P(I,J)) LINE3640 DU(I,J)=DU(I,J)/AP(I,J) LINE3650 103 CONTINUE LINE3660 CALL SOLVE LINE3670 100 CONTINUE LINE3680 CCOEFFICIENTS FOR THE V EQUATION------------------ LINE3690 NF=2 LINE3700 IF(.NOT.LSOLVE(NF)) GO TO 200 LINE3710 IST=2 LINE3720 JST=3 LINE3730 CALL GAMSOR LINE3740 REL=1.-RELAX(NF) LINE3750 DO 202 I=2,L2 LINE3760 AREA=R(1)*XCV(I) LINE3770 FLOW=AREA*V(I,2)*RHO(I,1) LINE3780 DIFF=AREA*GAM(I,1)/YCV(2) LINE3790 CALL DIFLOW LINE3800 202 AJM(I,3)=ACOF+AMAX1(0.,FLOW) LINE3810 DO 203 J=3,M2 LINE3820 FL=ARXJ(J)*U(2,J)*RHO(1,J) LINE3830 FLM=ARXJP(J-1)*U(2,J-1)*RHO(1,J-1) LINE3840 FLOW=FL+FLM LINE3850 DIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J)) LINE3860 CALL DIFLOW LINE3870 AIM(2,J)=ACOF+AMAX1(0.,FLOW) LINE3880 DO 203 I=2,L2 LINE3890 IF(I.EQ.L2) GO TO 204 LINE3900 FL=ARXJ(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) LINE3910 FLM=ARXJP(J-1)*U(I+1,J-1)*(FX(I+1)*RHO(I+1,J-1)+FXM(I+1)* LINE3920 1 RHO(I,J-1)) LINE3930 GM=GAM(I,J)*GAM(I+1,J)/(XCV(I)*GAM(I+1,J)+XCV(I+1)*GAM(I,J)+ LINE3940 1 1.E-30)*ARXJ(J) LINE3950 GMM=GAM(I,J-1)*GAM(I+1,J-1)/(XCV(I)*GAM(I+1,J-1)+XCV(I+1)* LINE3960 1 GAM(I,J-1)+1.0E-30)*ARXJP(J-1) LINE3970 DIFF=2.*(GM+GMM)/SXMN(J) LINE3980 GO TO 205 LINE3990 204 FL=ARXJ(J)*U(L1,J)*RHO(L1,J) LINE4000 FLM=ARXJP(J-1)*U(L1,J-1)*RHO(L1,J-1) LINE4010 DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J))LINE4020 205 FLOW=FL+FLM LINE4030 CALL DIFLOW LINE4040 AIM(I+1,J)=ACOF+AMAX1(0.,FLOW) LINE4050 AIP(I,J)=AIM(I+1,J)-FLOW LINE4060 IF(J.EQ.M2) GO TO 206 LINE4070 AREA=R(J)*XCV(I) LINE4080 FL=V(I,J)*(FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1))*RMN(J) LINE4090 FLP=V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))*RMN(J+1) LINE4100 FLOW=(FV(J)*FL+FVP(J)*FLP)*XCV(I) LINE4110 DIFF=AREA*GAM(I,J)/YCV(J) LINE4120 GO TO 207 LINE4130 206 AREA=R(M1)*XCV(I) LINE4140 FLOW=AREA*V(I,M1)*RHO(I,M1) LINE4150 DIFF=AREA*GAM(I,M1)/YCV(M2) LINE4160 207 CALL DIFLOW LINE4170 AJM(I,J+1)=ACOF+AMAX1(0.,FLOW) LINE4180 AJP(I,J)=AJM(I,J+1)-FLOW LINE4190 VOL=YCVRS(J)*XCV(I) LINE4200 SXT=SX(J) LINE4210 IF(J.EQ.M2) SXT=SX(M1) LINE4220 SXB=SX(J-1) LINE4230 IF(J.EQ.3) SXB=SX(1) LINE4240 APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)*RHO(I,J-1)* LINE4250 10.5*(SXB+SXMN(J)))/(YCVRS(J)*DT) LINE4260 AP(I,J)=AP(I,J)-APT LINE4270 CON(I,J)=CON(I,J)+APT*V(I,J) LINE4280 AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) LINE4290 1/RELAX(NF) LINE4300 CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*V(I,J) LINE4310 DV(I,J)=VOL/YDIF(J) LINE4320 CON(I,J)=CON(I,J)+DV(I,J)*(P(I,J-1)-P(I,J)) LINE4330 DV(I,J)=DV(I,J)/AP(I,J) LINE4340 203 CONTINUE LINE4350 CALL SOLVE LINE4360 200 CONTINUE LINE4370 CCOEFFICIENTS FOR THE PRESSURE CORRECTION EQUATION--------------------- LINE4380 NF=3 LINE4390 IF(.NOT.LSOLVE(NF)) GO TO 500 LINE4400 IST=2 LINE4410 JST=2 LINE4420 CALL GAMSOR LINE4430 SMAX=0. LINE4440 SSUM=0. LINE4450 DO 410 J=2,M2 LINE4460 DO 410 I=2,L2 LINE4470 VOL=YCVR(J)*XCV(I) LINE4480 410 CON(I,J)=CON(I,J)*VOL LINE4490 DO 402 I=2,L2 LINE4500 ARHO=R(1)*XCV(I)*RHO(I,1) LINE4510 CON(I,2)=CON(I,2)+ARHO*V(I,2) LINE4520 402 AJM(I,2)=0. LINE4530 DO 403 J=2,M2 LINE4540 ARHO=ARX(J)*RHO(1,J) LINE4550 CON(2,J)=CON(2,J)+ARHO*U(2,J) LINE4560 AIM(2,J)=0. LINE4570 DO 403 I=2,L2 LINE4580 IF(I.EQ.L2) GO TO 404 LINE4590 ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) LINE4600 FLOW=ARHO*U(I+1,J) LINE4610 CON(I,J)=CON(I,J)-FLOW LINE4620 CON(I+1,J)=CON(I+1,J)+FLOW LINE4630 AIP(I,J)=ARHO*DU(I+1,J) LINE4640 AIM(I+1,J)=AIP(I,J) LINE4650 GO TO 405 LINE4660 404 ARHO=ARX(J)*RHO(L1,J) LINE4670 CON(I,J)=CON(I,J)-ARHO*U(L1,J) LINE4680 AIP(I,J)=0. LINE4690 405 IF(J.EQ.M2) GO TO 406 LINE4700 ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J)) LINE4710 FLOW=ARHO*V(I,J+1) LINE4720 CON(I,J)=CON(I,J)-FLOW LINE4730 CON(I,J+1)=CON(I,J+1)+FLOW LINE4740 AJP(I,J)=ARHO*DV(I,J+1) LINE4750 AJM(I,J+1)=AJP(I,J) LINE4760 GO TO 407 LINE4770 406 ARHO=RMN(M1)*XCV(I)*RHO(I,M1) LINE4780 CON(I,J)=CON(I,J)-ARHO*V(I,M1) LINE4790 AJP(I,J)=0. LINE4800 407 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J) LINE4810 PC(I,J)=0. LINE4820 SMAX=AMAX1(SMAX,ABS(CON(I,J))) LINE4830 SSUM=SSUM+CON(I,J) LINE4840 403 CONTINUE LINE4850 CALL SOLVE LINE4860 CCOME HERE TO CORRECT THE PRESSURE AND VELOCITIES--------------------- LINE4870 DO 501 J=2,M2 LINE4880 DO 501 I=2,L2 LINE4890 P(I,J)=P(I,J)+PC(I,J)*RELAX(NP) LINE4900 IF(I.NE.2) U(I,J)=U(I,J)+DU(I,J)*(PC(I-1,J)-PC(I,J)) LINE4910 IF(J.NE.2) V(I,J)=V(I,J)+DV(I,J)*(PC(I,J-1)-PC(I,J)) LINE4920 501 CONTINUE LINE4930 500 CONTINUE LINE4940 CCOEFFICIENTS FOR OTHER EQUATIONS-------------------------------- LINE4950 IST=2 LINE4960 JST=2 LINE4970 DO 600 N=4,NFMAX LINE4980 NF=N LINE4990 IF(.NOT.LSOLVE(NF)) GO TO 600 LINE5000 CALL GAMSOR LINE5010 REL=1.-RELAX(NF) LINE5020 DO 602 I=2,L2 LINE5030 AREA=R(1)*XCV(I) LINE5040 FLOW=AREA*V(I,2)*RHO(I,1) LINE5050 DIFF=AREA*GAM(I,1)/YDIF(2) LINE5060 CALL DIFLOW LINE5070 602 AJM(I,2)=ACOF+AMAX1(0.,FLOW) LINE5080 DO 603 J=2,M2 LINE5090 FLOW=ARX(J)*U(2,J)*RHO(1,J) LINE5100 DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J)) LINE5110 CALL DIFLOW LINE5120 AIM(2,J)=ACOF+AMAX1(0.,FLOW) LINE5130 DO 603 I=2,L2 LINE5140 IF(I.EQ.L2) GO TO 604 LINE5150 FLOW=ARX(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) LINE5160 DIFF=ARX(J)*2.*GAM(I,J)*GAM(I+1,J)/((XCV(I)*GAM(I+1,J)+ LINE5170 1 XCV(I+1)*GAM(I,J)+1.0E-30)*SX(J)) LINE5180 GO TO 605 LINE5190 604 FLOW=ARX(J)*U(L1,J)*RHO(L1,J) LINE5200 DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J)) LINE5210 605 CALL DIFLOW LINE5220 AIM(I+1,J)=ACOF+AMAX1(0.,FLOW) LINE5230 AIP(I,J)=AIM(I+1,J)-FLOW LINE5240 AREA=RMN(J+1)*XCV(I) LINE5250 IF(J.EQ.M2) GO TO 606 LINE5260 FLOW=AREA*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J)) LINE5270 DIFF=AREA*2.*GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+ LINE5280 1 YCV(J+1)*GAM(I,J)+1.0E-30) LINE5290 GO TO 607 LINE5300 606 FLOW=AREA*V(I,M1)*RHO(I,M1) LINE5310 DIFF=AREA*GAM(I,M1)/YDIF(M1) LINE5320 607 CALL DIFLOW LINE5330 AJM(I,J+1)=ACOF+AMAX1(0.,FLOW) LINE5340 AJP(I,J)=AJM(I,J+1)-FLOW LINE5350 VOL=YCVR(J)*XCV(I) LINE5360 APT=RHO(I,J)/DT LINE5370 AP(I,J)=AP(I,J)-APT LINE5380 CON(I,J)=CON(I,J)+APT*F(I,J,NF) LINE5390 AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) LINE5400 1/RELAX(NF) LINE5410 CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*F(I,J,NF) LINE5420 603 CONTINUE LINE5430 CALL SOLVE LINE5440 600 CONTINUE LINE5450 TIME=TIME+DT LINE5460 ITER=ITER+1 LINE5470 IF(ITER.GE.LAST) LSTOP=.TRUE. LINE5480 RETURN LINE5490 END LINE5500 SUBROUTINE SUPPLY LINE5510 CHARACTER*10 TITLE LINE5520 LOGICAL LSOLVE,LPRINT,LBLK !,LSTOP LINE5530 PARAMETER (NOD=500) !NUMERO MAXIMO DE NODOS COMMON F(NOD,NOD,10),P(NOD,NOD),RHO(NOD,NOD),GAM(NOD,NOD), LINE5535 1 CON(NOD,NOD), LINE5540 2 AIP(NOD,NOD),AIM(NOD,NOD),AJP(NOD,NOD),AJM(NOD,NOD),AP(NOD,NOD), LINE5550 3 X(NOD),XU(NOD),XDIF(NOD),XCV(NOD),XCVS(NOD), LINE5560 4 Y(NOD),YV(NOD),YDIF(NOD),YCV(NOD),YCVS(NOD), LINE5570 5 YCVR(NOD),YCVRS(NOD),ARX(NOD),ARXJ(NOD),ARXJP(NOD), LINE5580 6 R(NOD),RMN(NOD),SX(NOD),SXMN(NOD),XCVI(NOD),XCVIP(NOD), LINE5590 7 YLS(10),M1S(10),M1ST(10),XLS(10),L1S(10),L1ST(10),NSX,NSY COMMON DU(NOD,NOD),DV(NOD,NOD),FV(NOD),FVP(NOD), LINE5600 1 FX(NOD),FXM(NOD),FY(NOD),FYM(NOD),PT(NOD),QT(NOD) LINE5610 COMMON/CHAR/TITLE(13) LINE5620 COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, LINE5630 1IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, LINE5640 2IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON LINE5650 REAL U(NOD,NOD),V(NOD,NOD),PC(NOD,NOD) LINE5660 C DIMENSION U(NOD,NOD),V(NOD,NOD),PC(NOD,NOD) LINE5660 EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1)) LINE5670 10 FORMAT(6(1H*),3X,A10,3X,6(1H*),/,8X,12(1H-)) LINE5680 20 FORMAT(4H I =,I6,12I9) LINE5690 30 FORMAT(1HJ) LINE5700 40 FORMAT(1X,I2,3X,(1x,13E9.2)) LINE5710 50 FORMAT(1H ) LINE5720 51 FORMAT(1X,'I =',2X,13(I4,5X)) LINE5730 52 FORMAT(1X,'X =',1x,13E9.2) LINE5740 53 FORMAT(1X,'TH =',1x,13E9.2) LINE5750 54 FORMAT(1X,'J =',2X,13(I4,5X)) LINE5760 55 FORMAT(1X,'Y =',1x,13E9.2) LINE5770 ENTRY UGRID LINE5780 XU(2)=0. LINE5790 DX=XL/FLOAT(L1-2) LINE5800 DO 1 I=3,L1 LINE5810 1 XU(I)=XU(I-1)+DX LINE5820 YV(2)=0. LINE5830 DY=YL/FLOAT(M1-2) LINE5840 DO 2 J=3,M1 LINE5850 2 YV(J)=YV(J-1)+DY LINE5860 RETURN LINE5870 ENTRY NUGRID L1ST(1)=L1S(1) DO 5 I=2,NSX 5 L1ST(I)=L1ST(I-1)+L1S(I) L1=L1ST(NSX) DO 6 I=1,NSX 6 XL=XLS(I)+XL XU(2)=0. I=1 DX=XLS(1)/FLOAT(L1S(1)-1) DO 3 J=3,L1 XU(J)=XU(J-1)+DX IF(J-1.EQ.L1ST(I)) THEN I=I+1 DX=XLS(I)/FLOAT(L1S(I)) IF(I.EQ.NSX) DX=XLS(I)/FLOAT(L1S(I)-1) 3 ENDIF M1ST(1)=M1S(1) DO 7 I=2,NSY 7 M1ST(I)=M1ST(I-1)+M1S(I) M1=M1ST(NSY) DO 8 I=1,NSY 8 YL=YLS(I)+YL YV(2)=0. I=1 DY=YLS(1)/FLOAT(M1S(1)-1) DO 4 J=3,M1 YV(J)=YV(J-1)+DY IF(J-1.EQ.M1ST(I)) THEN I=I+1 DY=YLS(I)/FLOAT(M1S(I)) IF(I.EQ.NSY) DY=YLS(I)/FLOAT(M1S(I)-1) 4 ENDIF RETURN ENTRY PRINT LINE5880 IF(.NOT.LPRINT(3)) GO TO 80 LINE5890 CALCULATE THE STREAM FUNCTION LINE5900 F(2,2,3)=0. LINE5910 DO 82 I=2,L1 LINE5920 IF(I.NE.2) F(I,2,3)=F(I-1,2,3)-RHO(I-1,1)*V(I-1,2) LINE5930 1*R(1)*XCV(I-1) LINE5940 DO 82 J=3,M1 LINE5950 RHOM=FX(I)*RHO(I,J-1)+FXM(I)*RHO(I-1,J-1) LINE5960 82 F(I,J,3)=F(I,J-1,3)+RHOM*U(I,J-1)*ARX(J-1) LINE5970 80 CONTINUE LINE5980 C LINE5990 IF(.NOT.LPRINT(NP)) GO TO 90 LINE6000 C LINE6010 CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATION LINE6020 DO 91 J=2,M2 LINE6030 P(1,J)=(P(2,J)*XCVS(3)-P(3,J)*XDIF(2))/XDIF(3) LINE6040 91 P(L1,J)=(P(L2,J)*XCVS(L2)-P(L3,J)*XDIF(L1))/XDIF(L2) LINE6050 DO 92 I=2,L2 LINE6060 P(I,1)=(P(I,2)*YCVS(3)-P(I,3)*YDIF(2))/YDIF(3) LINE6070 92 P(I,M1)=(P(I,M2)*YCVS(M2)-P(I,M3)*YDIF(M1))/YDIF(M2) LINE6080 P(1,1)=P(2,1)+P(1,2)-P(2,2) LINE6090 P(L1,1)=P(L2,1)+P(L1,2)-P(L2,2) LINE6100 P(1,M1)=P(2,M1)+P(1,M2)-P(2,M2) LINE6110 P(L1,M1)=P(L2,M1)+P(L1,M2)-P(L2,M2) LINE6120 PREF=P(IPREF,JPREF) LINE6130 DO 93 J=1,M1 LINE6140 DO 93 I=1,L1 LINE6150 93 P(I,J)=P(I,J)-PREF LINE6160 90 CONTINUE LINE6170 C LINE6180 c PRINT 50 LINE6190 c write(1,50) IEND=0 LINE6200 301 IF(IEND.EQ.L1) GO TO 310 LINE6210 IBEG=IEND+1 LINE6220 IEND=IEND+13 LINE6230 IEND=MIN0(IEND,L1) LINE6240 c PRINT 50 LINE6250 c write(1,50) c PRINT 51,(I,I=IBEG,IEND) LINE6260 c write(1,51)(i,i=ibeg,iend) IF(MODE.EQ.3)GO TO 302 LINE6270 c PRINT 52,(X(I),I=IBEG,IEND) LINE6280 c write(1,52)(X(I),I=IBEG,IEND) GO TO 303 LINE6290 302 continue c PRINT 53,(X(I),I=IBEG,IEND) LINE6300 c write(1,53)(X(I),I=IBEG,IEND) 303 GO TO 301 LINE6310 310 JEND=0 LINE6320 c PRINT 50 LINE6330 c write(1,50) 311 IF(JEND.EQ.M1) GO TO 320 LINE6340 JBEG=JEND+1 LINE6350 JEND=JEND+13 LINE6360 JEND=MIN0(JEND,M1) LINE6370 c PRINT 50 LINE6380 c write(1,50) c PRINT 54,(J,J=JBEG,JEND) LINE6390 c write(1,54)(J,J=JBEG,JEND) c PRINT 55,(Y(J),J=JBEG,JEND) LINE6400 c write(1,55)(y(j),J=JBEG,JEND) GO TO 311 LINE6410 320 CONTINUE LINE6420 C LINE6430 DO 999 N=1,NGAM LINE6440 NF=N LINE6450 IF(.NOT.LPRINT(NF))GO TO 999 LINE6460 c PRINT 50 LINE6470 c write(1,50) c PRINT 10,TITLE(NF) LINE6480 c write(1,10)title(nf) IFST=1 LINE6490 JFST=1 LINE6500 IF(NF.EQ.1.OR.NF.EQ.3) IFST=2 LINE6510 IF(NF.EQ.2.OR.NF.EQ.3) JFST=2 LINE6520 IBEG=IFST-13 LINE6530 110 CONTINUE LINE6540 IBEG=IBEG+13 LINE6550 IEND=IBEG+12 LINE6560 IEND=MIN0(IEND,L1) LINE6570 c PRINT 50 LINE6580 c write(1,50) c PRINT 20,(I,I=IBEG,IEND) LINE6590 c write(1,20)(I,I=IBEG,IEND) c PRINT 30 LINE6600 c write(1,30) JFL=JFST+M1 LINE6610 DO 115 JJ=JFST,M1 LINE6620 J=JFL-JJ LINE6630 c PRINT 40,J,(F(I,J,NF),I=IBEG,IEND) LINE6640 c write(1,40)J,(F(I,J,NF),I=IBEG,IEND) 115 CONTINUE LINE6650 IF(IEND.LT.L1) GO TO 110 LINE6660 999 CONTINUE LINE6670 RETURN LINE6680 END LINE6690 C LINE6700 C LINE6710