El modulo principal del programa de litodinamica para la zona costera c c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 c c Las ultimas correcciones: 18.07.01 C Mejoras: 12.02.02 c c--------------------------------------------------------------------------- INCLUDE 'COMML.FOR' INCLUDE 'VARS.INC' OPEN (1,FILE=LITOP//'ENLACE.DAT') READ(1,*)IEN CLOSE (1) OPEN (1,FILE=LITOP//'IPSI.DAT') READ(1,*)IPS READ(1,*)ITER1 READ(1,*)PERC CLOSE (1) OPEN (1,FILE=LITOP//'FORM.DAT') READ(1,*)IFORM READ(1,*)IFORMW READ(1,*)IFORMT CLOSE (1) IF(IEN.EQ.0)THEN C CALL INIC WRITE(*,*)' MODELO DE LITODINAMICA DE LA ZONA COSTERA (LIZC)' WRITE(*,*)' ' WRITE(*,*)'Centro de Investigaciones Oceanograficas e Hidrografica &s (CIOH)' WRITE(*,*)' ' WRITE(*,*)' Version 2.2, 1999 - 2001 ' WRITE(*,*)' ' WRITE(*,*)' ' ELSE WRITE(*,*)' ' WRITE(*,*)' ' WRITE(*,*)'MODELO DE LITODINAMICA' WRITE(*,*)' ' WRITE(*,*)'Continuaci�n. Enlace # ',IEN END IF CALL SWITCHE C------------------- GRID CONSTANTS --------------- OPEN (1,FILE=LITOP//'PARAM1.DAT') READ(1,*)HAG READ(1,*)PF READ(1,*)FI READ(1,*)M READ(1,*)N READ(1,*)MM READ(1,*)NN READ(1,*)DMIN READ(1,*)ICURV READ(1,*)SCALE READ(1,*)IXYP READ(1,*)IBOUND READ(1,*)XMIN READ(1,*)YMIN READ(1,*)XLIM READ(1,*)YLIM READ(1,*)MC,NC READ(1,*)NHC IF(NHC.GT.10)STOP'NO SE PERMITE LA CANTIDAD DE NODOS H(t) > 10' DO K=1,NHC READ(1,*)MCH(K),NCH(K) END DO C IF(ICURV.NE.0)HAG=1. ! puede ser en un caso especial YL=2.*7.29E-05*SIN(FI*3.14/180.) M1=M+1 N1=N+1 CLOSE (1) OPEN (1,FILE=LITOP//'PARAM3.DAT') READ(1,*)UNI IF(UNI.NE.0.)THEN READ(1,*)PF1 READ(1,*)PF2 READ(1,*)MREF READ(1,*)NREF OPEN (2,FILE=LITOP//'INDEX1.DAT') DO J=1,N1 IF(IFORM.EQ.1)JJ=N1-J+1 IF(IFORM.EQ.0)JJ=J READ(2,*)(ZNAK1(JJ,I),I=1,M1) END DO CLOSE (2) ELSE END IF CLOSE (1) C------------------- STRUCTURES ------------ OPEN (1,FILE=LITOP//'PARAM4.DAT') READ(1,*)PF3 CLOSE (1) OPEN (2,FILE=LITOP//'STRUCT.DAT') DO J=1,N1 IF(IFORM.EQ.1)JJ=N1-J+1 IF(IFORM.EQ.0)JJ=J READ(2,*)(ZNAK2(JJ,I),I=1,M1) END DO CLOSE (2) C------------------- DEPTH CALC. ------------------ OPEN (1,FILE=LITOP//'BOTTOM.DAT') DO J=1,N1 IF(IFORM.EQ.1)JJ=N1-J+1 IF(IFORM.EQ.0)JJ=J READ(1,*)(H(JJ,I),I=1,M1) END DO CLOSE (1) OPEN (1,FILE=LITOP//'BOTTOM0.DAT') ! DO J=1,N1 ! changeable name IF(IFORM.EQ.1)JJ=N1-J+1 ! IF(IFORM.EQ.0)JJ=J ! READ(1,*)(H0(JJ,I),I=1,M1) ! END DO ! CLOSE (1) ! CALL CHECK CALL BOT C----------------- PROFILES --------------------------- OPEN (1,FILE=LITOP//'PERFIL.DAT') READ(1,*)NPERF IF(NPERF.GT.10)THEN WRITE(*,*)'Cantidad de perfiles no puede ser mayor de 10' STOP ELSE END IF DO II=1,NPERF READ(1,*)MX(II),NY(II) IF(MX(II).NE.0.AND.NY(II).NE.0)THEN WRITE(*,*)'Coordenadas de los perfiles no son correctas' STOP ELSE END IF END DO READ(1,*)NPERF1 IF(NPERF1.GT.90)THEN WRITE(*,*)'Cantidad de perfiles digital no puede ser mayor de 90' STOP ELSE END IF IF(NPERF1.EQ.0)THEN CLOSE(1) ELSE READ(1,*)IAUTO READ(1,*)RMIN0 READ(1,*)RMIN1 READ(1,*)L5 DO II=1,NPERF1 READ(1,*)NP(II) DO JJ=1,NP(II) READ(1,*)XPP(II,JJ),YPP(II,JJ) END DO END DO CLOSE (1) END IF C------------------ INITIAL CONDS. --------------------- CALL ZEROS CALL PSINIT C--------------- METRIC TENSOR COMPONENTS ------------- IF(ICURV.EQ.0)GO TO 5 OPEN (9,FILE=RUTAP//'GRILLA/METRXY.OUT') READ(9,*)((x(J,I),i=1,m),j=1,n) READ(9,*)((y(J,I),i=1,m),j=1,n) READ(9,*)((x0(J,I),i=1,m1),j=1,n1) READ(9,*)((y0(J,I),i=1,m1),j=1,n1) CLOSE (9) c------------------------------------------------------------------ OPEN (8,FILE=RUTAP//'GRILLA/METRI.OUT') read(8,*)((HL1(J,I),I=1,M),J=1,N) read(8,*)((HL2(J,I),I=1,M),J=1,N) read(8,*)((HH1(J,I),I=1,M1),J=1,N1) read(8,*)((HH2(J,I),I=1,M1),J=1,N1) CLOSE (8) 5 CONTINUE C-------------- COMMON AND LITODYNAMIC PROPERTIES ---------------- OPEN(1,FILE=LITOP//'PARAM.DAT') READ(1,*)TYPE1 READ(1,*)IDAY READ(1,*)TTT READ(1,*)D50 READ(1,*)D90 READ(1,*)ROS READ(1,*)RO READ(1,*)Z0 READ(1,*)POROZ READ(1,*)ANU READ(1,*)DRIPLE READ(1,*)ICRIT READ(1,*)LEYF READ(1,*)UMOLIN READ(1,*)IORB READ(1,*)UORBM READ(1,*)HACTIV READ(1,*)TANTETA TANTETA=TANTETA*3.1415/180. GRAV=9.81 CLOSE (1) IF(ROS.LE.RO)STOP'SEDIMENT DENSITY IS LESS THAN THE WATER ONE' DCONST=((ROS/RO-1.)*GRAV/ANU/ANU)**(1./3.)*D50 DS=0.8*D50 ! SUSPENDED PARTICLE SIZE KSW=AMAX1(0.01,2.*DRIPLE) ! WAVE-RELATED RIPPLE BED ROUGHNESS Z0=AMAX1(0.01,Z0) ! CURRENT-RELATED BED ROUGHNESS CALL FALL CALL SHIELD c ADEAN=2.25*(WG*WG/GRAV)**(1./3.) ADEAN=0.067*(abs(WG)*100.)**0.44 c---------------------------------RESIDUAL TENSIONS --------------- IF(IRTF.EQ.0)GO TO 14 OPEN(1,FILE=LITOP//'RESIDS.DAT') DO J=1,N IF(IFORMT.EQ.1)JJ=N-J+1 IF(IFORMT.EQ.0)JJ=J READ(1,*)(FN1(JJ,I),I=1,M) END DO DO J=1,N IF(IFORMT.EQ.1)JJ=N-J+1 IF(IFORMT.EQ.0)JJ=J READ(1,*)(FN2(JJ,I),I=1,M) END DO DO J=1,N IF(IFORMT.EQ.1)JJ=N-J+1 IF(IFORMT.EQ.0)JJ=J READ(1,*)(U1NU(JJ,I),I=1,M) END DO DO J=1,N IF(IFORMT.EQ.1)JJ=N-J+1 IF(IFORMT.EQ.0)JJ=J READ(1,*)(U2NU(JJ,I),I=1,M) END DO CLOSE (1) 14 CONTINUE c------------------------------------------------------------------ OPEN (1,FILE=LITOP//'IMPRMR.DAT') ! READ THE KEYS FOR PRINT DO I=1,14 ! READ(1,*)IPRT(I) ! END DO ! CLOSE (1) c----------------------------------- SWAN RESULTS ----------------- IF(IWIC.EQ.0)GO TO 12 OPEN(1,FILE=SWANP//'FORCES.DAT') OPEN(2,FILE=SWANP//'HEIGHT.DAT') OPEN(3,FILE=SWANP//'LAMB.DAT') OPEN(4,FILE=SWANP//'TAU.DAT') OPEN(5,FILE=SWANP//'DIRECT.DAT') DO J=1,N IF(IFORMW.EQ.1)JJ=N-J+1 IF(IFORMW.EQ.0)JJ=J READ(1,*)(FX(JJ,I),I=1,M) READ(2,*)(HW(JJ,I),I=1,M) READ(3,*)(FLAMB(JJ,I),I=1,M) READ(4,*)(TAU(JJ,I),I=1,M) READ(5,*)(DIR(JJ,I),I=1,M) END DO DO J=1,N IF(IFORMW.EQ.1)JJ=N-J+1 IF(IFORMW.EQ.0)JJ=J READ(1,*)(FY(JJ,I),I=1,M) END DO CLOSE (1) CLOSE (2) CLOSE (3) CLOSE (4) CLOSE (5) OPEN(1,FILE=SWANP//'ORBIT.DAT') DO J=1,N IF(IFORMW.EQ.1)JJ=N-J+1 IF(IFORMW.EQ.0)JJ=J READ(1,*)(UORB(JJ,I),I=1,M) END DO CLOSE (1) CALL TESTWA IF(ISMOOTW.EQ.1) CALL SMOOTHW 12 CONTINUE c------------------------------------------------------------------ OPEN(10,FILE=LITOP//'WIND.DAT') READ(10,*)WINDX,WINDY CALL WINDWAV ! FOR CONSTANT WINDS ONLY (SEE RHS) CLOSE (10) C------------------------------------------------------ IF(MM.GT.M)STOP'El| area reducida es mayor que el area completa' IF(NN.GT.N)STOP'El area reducida es mayor que el area completa' IF(IEN.EQ.0)THEN IF(MM.LT.M.OR.NN.LT.N)WRITE(*,*)'ATENCION! El area de calculo esta| & reducida' ELSE END IF M0=M N0=N M=MM N=NN M1=M+1 N1=N+1 C------------------------------------------------------ INICIO ..... CALL PRNH CALL SLOPE IF(TYPE1.EQ.1)IDAY=1 WRITE(*,*)' ' WRITE(*,*)' ' IF(IEN.EQ.0)THEN WRITE(*,*)'INICIO DE CALCULO ...' WRITE(*,*)' ' WRITE(*,*)' ' IF(IPS.EQ.1)THEN OPEN(1,FILE=LITOP//'PSI0.DAT') READ(1,*)((PSIN(J,I),J=1,N1),I=1,M1) DO J=1,N1 DO I=1,M1 PSI(J,I)=PSIN(J,I) PSI0(J,I)=PSIN(J,I) END DO END DO CLOSE (1) ELSE END IF ELSE END IF DO 100 KT=1,IDAY ! 12 MESES DE CALCULO CALL VORTX1 WRITE(*,33)KT,L 33 FORMAT('+CANTIDAD DE PASOS Y DE ITERACIONES = ',2I8) CALL WAVEL CALL PROFU CALL PROFK CALL PROFC IF(IWIC.NE.0)CALL WAVEUP IF(ICRIT.EQ.84)THEN CALL UCRIT84 CALL FLUX84 ELSE CALL UCRIT90 CALL FLUX90 END IF CALL FONDO CALL COSTA ! CALL OUTPUT 100 CONTINUE IF(IPS.EQ.1)THEN OPEN(1,FILE=LITOP//'PSI0.DAT') WRITE(1,*)((PSI(J,I),J=1,N1),I=1,M1) CLOSE (1) ELSE END IF C-------------------------------- OUTPUT ---------------- C CALL OUTPUT C CALL FINAL IF(TYPE1.EQ.2)STOP WRITE(*,*)' ' write(*,*)'El calculo fue terminado exitosamente' pause STOP END SUBROUTINE VORTX1 c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' epsp=0.0001 MOL=0 DO 19 I=1,M1 DO 19 J=1,N1 IF(KT.EQ.1)PSI(J,I)=PSI0(J,I) IF(H(J,I).EQ.PF)GOTO 19 MOL=MOL+1 19 CONTINUE TETA=1. L=0 35 L=L+1 IF(L.GE.2000000)GOTO 325 CALL BOUND DO 778 I=2,M DO 778 J=2,N IF(H(J,I).EQ.PF)GOTO 778 C IF(H(J+1,I).EQ.PF.OR.H(J,I+1).EQ.PF.OR.H(J,I-1).EQ.PF C &.OR.H(J-1,I).EQ.PF)GOTO 778 C IF(H(J+1,I+1).EQ.PF.OR.H(J-1,I+1).EQ.PF.OR.H(J+1,I-1).EQ.PF C &.OR.H(J-1,I-1).EQ.PF)GOTO 778 DDT=(AFV(J,I)+AFV(J-1,I)-AFV(J,I-1)-AFV(J-1,I-1) &-AFU(J,I)-AFU(J,I-1)+AFU(J-1,I)+AFU(J-1,I-1)) A1=RK(J,I)*HL1(J,I)/HL2(J,I)/H1(J,I) A4=RK(J,I-1)*HL1(J,I-1)/HL2(J,I-1)/H1(J,I-1) A8=RK(J-1,I)*HL1(J-1,I)/HL2(J-1,I)/H1(J-1,I) A9=RK(J-1,I-1)*HL1(J-1,I-1)/HL2(J-1,I-1)/H1(J-1,I-1) B4=YL*HH1(J+1,I)/HH2(J+1,I)/H(J+1,I) B8=YL*HH1(J-1,I)/HH2(J-1,I)/H(J-1,I) G1=RK(J,I)*HL2(J,I)/HL1(J,I)/H1(J,I) G2=RK(J-1,I)*HL2(J-1,I)/HL1(J-1,I)/H1(J-1,I) G6=RK(J,I-1)*HL2(J,I-1)/HL1(J,I-1)/H1(J,I-1) G7=RK(J-1,I-1)*HL2(J-1,I-1)/HL1(J-1,I-1)/H1(J-1,I-1) O2=YL*HH2(J,I+1)/HH1(J,I+1)/H(J,I+1) O6=YL*HH2(J,I-1)/HH1(J,I-1)/H(J,I-1) PSI1=((G2+G1)*PSI0(J,I+1)+(G7+G6)*PSI(J,I-1)+ & (A4+A1)*PSI0(J+1,I)+(A9+A8)*PSI(J-1,I)+ & 0.5*(B4*(PSI0(J+1,I+1)-PSI0(J+1,I-1))-B8*(PSI0(J-1,I+1)- & PSI(J-1,I-1))-O2*(PSI0(J+1,I+1)-PSI0(J-1,I+1))+ & O6*(PSI0(J+1,I-1)-PSI(J-1,I-1)))-HAG*DDT)/ ! DOESN'T DIVIDE BY 2 AND ONE HAG! & (G1+G2+G6+G7+A1+A4+A8+A9) PSI(J,I)=TETA*PSI1+(1.-TETA)*PSI0(J,I) 778 CONTINUE CALL VELCTY CALL FRIC CALL RHS if(L.eq.1)go to 35 KAP=0 DO 97 I=1,M1 DO 97 J=1,N1 IF(H(J,I).EQ.PF)GOTO 97 IF(PSI0(J,I).EQ.0..AND.PSI(J,I).EQ.0.)GOTO 322 IF(ABS((PSI(J,I)-PSI0(J,I))/(PSI0(J,I)+0.1E-09)).GT.EPSP)GOTO 97 322 KAP=KAP+1 97 CONTINUE DO 2 I=1,M1 DO 2 J=1,N1 IF(H(J,I).EQ.PF)GOTO 2 PSI0(J,I)=PSI(J,I) 2 CONTINUE IF(KAP.EQ.MOL)GOTO 444 IF(L/ITER1*ITER1.EQ.L)WRITE(*,*) &'Acercamiento esta en ',KAP/FLOAT(MOL)*100,' %' IF(KAP/FLOAT(MOL)*100..GE.PERC)GO TO 444 GOTO 35 444 CONTINUE CCC write(*,*)l,kap,mol RETURN 325 WRITE(*,764) 764 FORMAT(2X,'THE NUMBER OF ITERATIONS FOR PSI IS MORE THEN 2e4'/) RETURN END SUBROUTINE BOUND0 c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' c---------------- eastern and western bounds ----- DO 1 J=2,N if(h(j,M1).eq.pf)go to 1 c PSI(J,M1)=PSI(j,m)*2.-psi(j,m-1) c PSI0(J,M1)=PSI(j,m)*2.-psi(j,m-1) PSI(J,M1)=PSI(j,m) PSI0(J,M1)=PSI(j,m) 1 continue DO 4 J=2,N if(h(j,1).eq.pf)go to 4 c PSI(J,1)=PSI(j,2)*2.-psi(j,3) c PSI0(J,1)=PSI(j,2)*2.-psi(j,3) PSI(J,1)=PSI(j,2) PSI0(J,1)=PSI(j,2) 4 continue c---------------- southern-nothern bounds ----- DO 3 i=2,M if(h(1,i).eq.pf)go to 2 c PSI(1,I)=PSI(2,I)*2.-PSI(3,I) c PSI0(1,I)=PSI(2,I)*2.-PSI(3,I) PSI(1,I)=PSI(2,I) PSI0(1,I)=PSI(2,I) 2 continue if(h(N1,i).eq.pf)go to 3 PSI(N1,I)=PSI(N,I) PSI0(N1,I)=PSI(N,I) 3 continue C--------------- COSMETICS -------------------- IF(H(1,M1).NE.PF)PSI(1,M1)=PSI(1,M)+PSI(2,M1)-PSI(2,M) IF(H(N1,M1).NE.PF)PSI(N1,M1)=PSI(N1,M)+PSI(N,M1)-PSI(N,M) IF(H(N1,1).NE.PF)PSI(N1,1)=PSI(N1,2)+PSI(N,1)-PSI(N,2) IF(H(1,1).NE.PF)PSI(1,1)=PSI(1,2)+PSI(2,1)-PSI(2,2) IF(H(1,M1).NE.PF)PSI0(1,M1)=PSI(1,M)+PSI(2,M1)-PSI(2,M) IF(H(N1,M1).NE.PF)PSI0(N1,M1)=PSI(N1,M)+PSI(N,M1)-PSI(N,M) IF(H(N1,1).NE.PF)PSI0(N1,1)=PSI(N1,2)+PSI(N,1)-PSI(N,2) IF(H(1,1).NE.PF)PSI0(1,1)=PSI(1,2)+PSI(2,1)-PSI(2,2) C------------------------------------------------- IF(UNI.EQ.0.)RETURN DO 10 I=1,M1 DO 10 J=1,N1 IF(ZNAK1(J,I).EQ.PF1)THEN PSI(J,I)=0. PSI0(J,I)=0. ELSE END IF 10 CONTINUE IF(MREF.EQ.0)GO TO 50 C A LO LARGO DEL EJE Y I=MREF DO 21 J=1,N1 IF(ZNAK1(J,I).EQ.PF1)NREF1=J 21 CONTINUE VALOR=0. DO 41 J=NREF1,N1 IF(H(J+1,I).EQ.PF)GO TO 41 VALOR=VALOR-U(J,I)*HAG*H1(J,I)*HL1(J,I)*HL2(J,I) 41 CONTINUE GO TO 60 C C A LO LARGO DEL EJE X 50 J=NREF DO 20 I=1,M1 IF(ZNAK1(J,I).EQ.PF1)MREF1=I 20 CONTINUE VALOR=0. DO 40 I=MREF1,M1 IF(H(J,I+1).EQ.PF)GO TO 40 VALOR=VALOR+V(J,I)*HAG*H1(J,I)*HL1(J,I)*HL2(J,I) 40 CONTINUE 60 CONTINUE C DO 30 I=1,M1 DO 30 J=1,N1 IF(ZNAK1(J,I).EQ.PF2)THEN PSI(J,I)=VALOR PSI0(J,I)=VALOR ELSE END IF 30 CONTINUE RETURN END SUBROUTINE VELCTY c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.0, 2000 INCLUDE 'COMML.FOR' DO 174 I=1,M DO 174 J=1,N u(j,i)=0. v(j,i)=0. IF(HL(J,I).EQ.PF)GOTO 174 QA=H1(J,I)*HL1(J,I)*HL2(J,I) U(J,I)=-(PSI(J+1,I+1)+PSI(J+1,I)-PSI(J,I+1)-PSI(J,I))/(2.*HAG*QA)- & U1NU(J,I)/H1(J,I) V(J,I)=(PSI(J+1,I+1)+PSI(J,I+1)-PSI(J+1,I)-PSI(J,I))/(2.*HAG*QA)- & U2NU(J,I)/H1(J,I) 174 CONTINUE RETURN END subroutine tang(b1,a1,c) p=3.14159265 if(b1.eq.0.)c=p/2. if(b1.ne.0.)c=abs(atan(a1/b1)) l=sign(1.,a1) l=l+1 n=sign(1.,b1) n=n+1 if(l.eq.2.and.n.eq.0)c=p-c if(l.eq.0.and.n.eq.0)c=p+c if(l.eq.0.and.n.eq.2)c=2*p-c C=C*180./P return end SUBROUTINE UCRIT84 c Elaborado por Serguei A. Lonin, 1999-2000 Version 1.1, 2000 INCLUDE 'COMML.FOR' IF(D50.LT.1.E-04)STOP'THERE ARE VERY SMALL PARTICLES' IF(D50.GT.2.E-03)STOP'THERE ARE VERY LARGE PARTICLES' IF(D50.LT.5.E-04)THEN DO 1 J=1,N DO 1 I=1,M UCR(J,I)=1.E10 IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 UCR(J,I)=0.19*D50**0.1*ALOG10(4.*H1(J,I)/D90) 1 CONTINUE ELSE DO 2 J=1,N DO 2 I=1,M UCR(J,I)=1.E10 IF(HL(J,I).EQ.PF)GO TO 2 IF(H1(J,I).LT.DMIN)GO TO 2 UCR(J,I)=8.5*D50**0.6*ALOG10(4.*H1(J,I)/D90) 2 CONTINUE END IF IF(IWIC.NE.0)CALL WAVECR RETURN END SUBROUTINE UCRIT90 c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.0, 2000 INCLUDE 'COMML.FOR' DO 1 J=1,N DO 1 I=1,M UCR(J,I)=1.E10 IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 UCR(J,I)=5.75*SQRT((ROS/RO-1.)*GRAV*D50*TETACR)* & ALOG10(4.*H1(J,I)/D90) 1 CONTINUE IF(IWIC.NE.0)CALL WAVECR RETURN END SUBROUTINE WAVECR c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' DO 1 J=1,N DO 1 I=1,M UCRW(J,I)=1.E10 ! ? SEE ZEROS ? IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 IF(TAU(J,I).LE.0.)GO TO 1 IF(D50.LT.0.0005)THEN UCRW(J,I)=(0.12*(ROS/RO-1.)*GRAV*SQRT(D50)*SQRT(TAU(J,I))) & **(2./3.) ELSE UCRW(J,I)=(1.09*(ROS/RO-1.)*GRAV*(D50)**0.75*TAU(J,I)**0.25) & **0.571 END IF 1 CONTINUE RETURN END SUBROUTINE FLUX84 c Elaborado por Serguei A. Lonin, 1999-2000 Version 1.1, 2000 INCLUDE 'COMML.FOR' DO 1 J=1,N DO 1 I=1,M QSX(J,I)=0. ! QSY(J,I)=0. ! NOT SO FOR ABRASION CASE QBX(J,I)=0. ! QBY(J,I)=0. ! IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 UEFF=SQRT(U(J,I)*U(J,I)+V(J,I)*V(J,I))+UORB(J,I) ! TRANS=(AMAX1(0.,(UEFF-UCR(J,I)))/SQRT((ROS/RO-1.)*GRAV*D50))**2.4 QSX(J,I)=0.012*U(J,I)*H1(J,I)*TRANS*D50/H1(J,I)/DCONST**0.6 & *ROS*IQS QSY(J,I)=0.012*V(J,I)*H1(J,I)*TRANS*D50/H1(J,I)/DCONST**0.6 & *ROS*IQS QBX(J,I)=0.005*U(J,I)*H1(J,I)*TRANS*(D50/H1(J,I))**1.2 & *ROS*IQB QBY(J,I)=0.005*V(J,I)*H1(J,I)*TRANS*(D50/H1(J,I))**1.2 & *ROS*IQB QTX(J,I)=QSX(J,I)+QBX(J,I) QTY(J,I)=QSY(J,I)+QBY(J,I) 1 CONTINUE RETURN END SUBROUTINE FONDO c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' ACONST=TTT/2./HAG/ROS/(1.-POROZ) CALL LIMIT_S DO 1 I=2,M DO 1 J=2,N DELTA(J,I)=0. IF(H(J,I).EQ.PF)GO TO 1 IF(H(J,I).LT.DMIN)GO TO 1 DIV=HL1(J,I)*QTX(J,I)-HL1(J,I-1)*QTX(J,I-1)+HL1(J-1,I)*QTX(J-1,I) & -HL1(J-1,I-1)*QTX(J-1,I-1) & +HL2(J,I)*QTY(J,I)-HL2(J-1,I)*QTY(J-1,I)+HL2(J,I-1)*QTY(J,I-1)- & HL2(J-1,I-1)*QTY(J-1,I-1) DELTA(J,I)=ACONST/HH1(J,I)/HH2(J,I)*DIV *0.1 ! lony IF(TYPE1.EQ.2)THEN DELTA(J,I)=AMAX1(DELTA(J,I),DLIMM) DELTA(J,I)=AMIN1(DELTA(J,I),DLIMP) HP(J,I)=H(J,I) H(J,I)=H(J,I)+DELTA(J,I) IF(H(J,I).LT.DMIN)THEN H(J,I)=PF ELSE END IF ELSE END IF 1 CONTINUE IF(TYPE1.EQ.2)THEN IF(ISMOOT.EQ.1)THEN DO ISM=1,KRATS CALL SMOOTH END DO DO 11 I=2,M DO 11 J=2,N IF(H(J,I).EQ.PF)GO TO 11 IF(H(J,I).LT.DMIN)GO TO 11 DELTA(J,I)=H(J,I)-HP(J,I) 11 CONTINUE ELSE END IF CALL BOUNDH ! was opened on September 2001 CALL BOT ELSE END IF RETURN END SUBROUTINE BOT c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' DO 107 I=1,M DO 107 J=1,N FT1=H(J,I) FT2=H(J+1,I) FT3=H(J,I+1) FT4=H(J+1,I+1) NUMH=4 IF(FT1.EQ.PF)THEN FT1=0. NUMH=NUMH-1 ELSE END IF IF(FT2.EQ.PF)THEN FT2=0. NUMH=NUMH-1 ELSE END IF IF(FT3.EQ.PF)THEN FT3=0. NUMH=NUMH-1 ELSE END IF IF(FT4.EQ.PF)THEN FT4=0. NUMH=NUMH-1 ELSE END IF IF(NUMH.NE.0)THEN H1(J,I)=(FT1+FT2+FT3+FT4)/FLOAT(NUMH) ELSE H1(J,I)=0. END IF C H1(J,I)=AMin1(ft1,ft2,ft3,ft4) HL(J,I)=AMin1(H(J,I),H(J+1,I),H(J,I+1),H(J+1,I+1)) ccccccHL(J,I)=AMAX1(H(J,I),H(J+1,I),H(J,I+1),H(J+1,I+1)) ! error !!?? 107 CONTINUE RETURN END SUBROUTINE BOUNDH c Elaborado por Serguei A. Lonin, 1999-2000 Version 1.1, 2000 INCLUDE 'COMML.FOR' C----------------- EASTERN AND WESTERN BOUNDS ----- do j=2,n IF(U(J,M).NE.0.)h(j,m1)=h(j,m) IF(U(J,1).NE.0.)H(J,1)=H(J,2) end do C----------------- NORTHERN AND SOUTHERN BOUNDS ----- do I=2,M IF(V(N,I).NE.0.)H(N1,I)=H(N,I) IF(V(1,I).NE.0.)H(1,I)=H(2,I) end do RETURN END SUBROUTINE FRIC c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' RAD=3.14/180. DO 174 I=1,M DO 174 J=1,N IF(HL(J,I).EQ.PF)GOTO 174 CD=1./32./(ALOG10(14.8*AMAX1(DMIN,H1(J,I))/Z0))**2 UMO=UMOLIN IF(LEYF.EQ.0)GO TO 1 UMO=SQRT(U(J,I)*U(J,I)+V(J,I)*V(J,I)) CALL TANG(U(J,I),V(J,I),ANGU) IF(IORB.EQ.0)THEN UMO=SQRT(UMO**2+UORBM**2) ELSE UMO=SQRT(UMO**2+UORB(J,I)**2+2.*UMO*UORB(J,I)*COS(RAD* & (ANGU-DIR(J,I)))) END IF 1 UMO=AMAX1(UMO,0.001) RK(J,I)=CD*UMO/AMAX1(DMIN,H1(J,I)) 174 CONTINUE RETURN END SUBROUTINE PSINIT c Elaborado por Serguei A. Lonin, 1999-2000 Version 1.1, 2000 C Este programa debe ser dise�ada por el usuario INCLUDE 'COMML.FOR' INCLUDE 'VARS.INC' DO I=1,M1 DO J=1,N1 PSI0(J,I)=0. END DO END DO OPEN(1,FILE=LITOP//'PSI0.DAT') WRITE(1,*)((PSI0(J,I),J=1,N1),I=1,M1) CLOSE (1) C DO 2 I=1,M1 C DO 1 J=1,N1 C IF(J.LT.40)THEN C PSI0(J,I)=0. C ELSE C PSI0(J,I)=14.54 ! CAUDAL C END IF C 1 CONTINUE C 2 CONTINUE C QTX(40,22)=-0.03 ! SEDIMENT FLUX C return end SUBROUTINE WINDWAV c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' C---------------- INPUT PHYSICAL WIND STRESS ------------- DO 50 I=1,M DO 50 J=1,N T1(J,I)=0. T2(J,I)=0. FM1(J,I)=0. FM2(J,I)=0. 50 CONTINUE WINDM=1.25*1.3E-03*SQRT(WINDX*WINDX+WINDY*WINDY) TAUX=WINDM*WINDX*IWC TAUY=WINDM*WINDY*IWC C------------------------------------------- DO 100 I=2,M-1 DO 100 J=2,N-1 YN=(Y(J+1,I)-Y(J-1,I))*0.5/HAG YK=(Y(J,I+1)-Y(J,I-1))*0.5/HAG XN=(X(J+1,I)-X(J-1,I))*0.5/HAG XK=(X(J,I+1)-X(J,I-1))*0.5/HAG AK=XK*YN-XN*YK IF(AK.EQ.0.)GO TO 100 T1(J,I)=(TAUX*YN-TAUY*XN)/AK T2(J,I)=(TAUY*XK-TAUX*YK)/AK FM1(J,I)=(FX(J,I)*YN-FY(J,I)*XN)/AK FM2(J,I)=(FY(J,I)*XK-FX(J,I)*YK)/AK 100 CONTINUE c------------------------------------------------------------------ DO 101 I=2,M-1 J=1 YN=(Y0(J+1,I)-Y0(J,I)+Y0(J+1,I+1)-Y0(J,I+1))*0.5/HAG YK=(Y0(J,I+1)-Y0(J,I)+Y0(J+1,I+1)-Y0(J+1,I))*0.5/HAG XN=(X0(J+1,I)-X0(J,I)+X0(J+1,I+1)-X0(J,I+1))*0.5/HAG XK=(X0(J,I+1)-X0(J,I)+X0(J+1,I+1)-X0(J+1,I))*0.5/HAG AK=XK*YN-XN*YK IF(AK.EQ.0.)GO TO 109 T1(J,I)=(TAUX*YN-TAUY*XN)/AK T2(J,I)=(TAUY*XK-TAUX*YK)/AK FM1(J,I)=(FX(J,I)*YN-FY(J,I)*XN)/AK FM2(J,I)=(FY(J,I)*XK-FX(J,I)*YK)/AK 109 J=N YN=(Y0(J+1,I)-Y0(J,I)+Y0(J+1,I+1)-Y0(J,I+1))*0.5/HAG YK=(Y0(J,I+1)-Y0(J,I)+Y0(J+1,I+1)-Y0(J+1,I))*0.5/HAG XN=(X0(J+1,I)-X0(J,I)+X0(J+1,I+1)-X0(J,I+1))*0.5/HAG XK=(X0(J,I+1)-X0(J,I)+X0(J+1,I+1)-X0(J+1,I))*0.5/HAG AK=XK*YN-XN*YK IF(AK.EQ.0.)GO TO 101 T1(J,I)=(TAUX*YN-TAUY*XN)/AK T2(J,I)=(TAUY*XK-TAUX*YK)/AK FM1(J,I)=(FX(J,I)*YN-FY(J,I)*XN)/AK FM2(J,I)=(FY(J,I)*XK-FX(J,I)*YK)/AK 101 CONTINUE DO 102 J=2,N-1 I=1 YN=(Y0(J+1,I)-Y0(J,I)+Y0(J+1,I+1)-Y0(J,I+1))*0.5/HAG YK=(Y0(J,I+1)-Y0(J,I)+Y0(J+1,I+1)-Y0(J+1,I))*0.5/HAG XN=(X0(J+1,I)-X0(J,I)+X0(J+1,I+1)-X0(J,I+1))*0.5/HAG XK=(X0(J,I+1)-X0(J,I)+X0(J+1,I+1)-X0(J+1,I))*0.5/HAG AK=XK*YN-XN*YK IF(AK.EQ.0.)GO TO 108 T1(J,I)=(TAUX*YN-TAUY*XN)/AK T2(J,I)=(TAUY*XK-TAUX*YK)/AK FM1(J,I)=(FX(J,I)*YN-FY(J,I)*XN)/AK FM2(J,I)=(FY(J,I)*XK-FX(J,I)*YK)/AK 108 I=M YN=(Y0(J+1,I)-Y0(J,I)+Y0(J+1,I+1)-Y0(J,I+1))*0.5/HAG YK=(Y0(J,I+1)-Y0(J,I)+Y0(J+1,I+1)-Y0(J+1,I))*0.5/HAG XN=(X0(J+1,I)-X0(J,I)+X0(J+1,I+1)-X0(J,I+1))*0.5/HAG XK=(X0(J,I+1)-X0(J,I)+X0(J+1,I+1)-X0(J+1,I))*0.5/HAG AK=XK*YN-XN*YK IF(AK.EQ.0.)GO TO 102 T1(J,I)=(TAUX*YN-TAUY*XN)/AK FM1(J,I)=(FX(J,I)*YN-FY(J,I)*XN)/AK FM2(J,I)=(FY(J,I)*XK-FX(J,I)*YK)/AK 102 CONTINUE T1(1,1)=T1(1,2)+T1(2,1)-T1(2,2) T1(N,1)=T1(N,2)+T1(N-1,1)-T1(N-1,2) T1(1,M)=T1(1,M-1)+T1(2,M)-T1(2,M-1) T1(N,M)=T1(N,M-1)+T1(N-1,M)-T1(N-1,M-1) T2(1,1)=T2(1,2)+T2(2,1)-T2(2,2) T2(N,1)=T2(N,2)+T2(N-1,1)-T2(N-1,2) T2(1,M)=T2(1,M-1)+T2(2,M)-T2(2,M-1) T2(N,M)=T2(N,M-1)+T2(N-1,M)-T2(N-1,M-1) FM1(1,1)=FM1(1,2)+FM1(2,1)-FM1(2,2) FM1(N,1)=FM1(N,2)+FM1(N-1,1)-FM1(N-1,2) FM1(1,M)=FM1(1,M-1)+FM1(2,M)-FM1(2,M-1) FM1(N,M)=FM1(N,M-1)+FM1(N-1,M)-FM1(N-1,M-1) FM2(1,1)=FM2(1,2)+FM2(2,1)-FM2(2,2) FM2(N,1)=FM2(N,2)+FM2(N-1,1)-FM2(N-1,2) FM2(1,M)=FM2(1,M-1)+FM2(2,M)-FM2(2,M-1) FM2(N,M)=FM2(N,M-1)+FM2(N-1,M)-FM2(N-1,M-1) DO 51 I=1,M DO 51 J=1,N T1(J,I)=T1(J,I)/RO T2(J,I)=T2(J,I)/RO IF(HL(J,I).EQ.PF)GO TO 51 FM1(J,I)=-FM1(J,I)/RO/H1(J,I) ! - FM2(J,I)=-FM2(J,I)/RO/H1(J,I) ! - 51 CONTINUE RETURN END SUBROUTINE RHS c Elaborado por Serguei A. Lonin, 1999-2000 Version 1.1, 2000 INCLUDE 'COMML.FOR' C FM1 & FM2 - AVERAGE WAVE ACTION IN RHS C FL1 & FL2 - AVERAGE NON LINEAR FORCES C FN1 & FN2 - AVERAGE NIHOUL'S TERMS (NO MODIFIED) C FN1MD & FN2MD - AVERAGE NIHOUL'S TERMS (MODIFIED) IF(INLT.NE.0)CALL NOLIN DO 1 I=1,M DO 1 J=1,N IF(HL(J,I).EQ.PF)GO TO 1 FN1MD(J,I)=FN1(J,I)-(YL*U2NU(J,I)+RK(J,I)*U1NU(J,I))/H1(J,I) FN2MD(J,I)=FN2(J,I)-(-YL*U1NU(J,I)+RK(J,I)*U2NU(J,I))/H1(J,I) RHSX=FN1MD(J,I)-FM1(J,I)-FL1(J,I) RHSY=FN2MD(J,I)-FM2(J,I)-FL2(J,I) AFU(J,I)=HL1(J,I)*HL1(J,I)*(T1(J,I)/H1(J,I)+RHSX) AFV(J,I)=HL2(J,I)*HL2(J,I)*(T2(J,I)/H1(J,I)+RHSY) 1 CONTINUE RETURN END SUBROUTINE NOLIN c Elaborado por Serguei A. Lonin, 1999-2000 Version 1.1, 2000 INCLUDE 'COMML.FOR' DO 23 I=2,M-1 DO 23 J=2,N-1 IF(HL(J,I).EQ.PF)GO TO 23 G1=HL1(J,I+1)*HL1(J,I+1) G2=HL1(J,I-1)*HL1(J,I-1) G3=HL1(J+1,I)*HL1(J+1,I) G5=HL2(J,I+1)*HL2(J,I+1) G6=HL2(J,I-1)*HL2(J,I-1) G7=HL2(J+1,I)*HL2(J+1,I) G01=HL1(J,I)*HL1(J,I) G02=HL2(J,I)*HL2(J,I) G4=HL1(J-1,I)*HL1(J-1,I) G8=HL2(J-1,I)*HL2(J-1,I) FL1(J,I)=((U(J,I)+U(J,I+1))*U(J,I+1)-(U(J,I-1)+ & U(J,I))*U(J,I-1)+(V(J+1,I)+V(J,I))*U(J+1,I)- & (V(J-1,I)+V(J,I))*U(J-1,I))/HAG*0.25+ & (U(J,I)*U(J,I)*0.5*(ALOG(G1)-ALOG(G2))+ & U(J,I)*V(J,I)*(ALOG(G3)-ALOG(G4))- & V(J,I)*V(J,I)*0.5/G01*(G5-G6))*0.5/HAG FL2(J,I)=((U(J,I)+U(J,I+1))*V(J,I+1)-(U(J,I-1)+ & U(J,I))*V(J,I-1)+(V(J+1,I)+V(J,I))*V(J+1,I)- & (V(J-1,I)+V(J,I))*V(J-1,I))/HAG*0.25- & (U(J,I)*U(J,I)*0.5/G02*(G3-G4)+ & U(J,I)*V(J,I)*(ALOG(G5)-ALOG(G6))+ & V(J,I)*V(J,I)*0.5*(ALOG(G7)-ALOG(G8)))*0.5/HAG 23 CONTINUE RETURN END SUBROUTINE WAVEUP c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' C FM1 & FM2 - AVERAGE WAVE ACTION IN RHS C FL1 & FL2 - AVERAGE NON LINEAR FORCES C FN1 & FN2 - AVERAGE NIHOUL'S TERMS (NO MODIFIED) DO 1 I=1,M DO 1 J=1,N IF(HL(J,I).EQ.PF)GO TO 1 GRAD1(J,I)=FN1(J,I)-FM1(J,I)+T1(J,I)/H1(J,I)-FL1(J,I)+YL*V(J,I)- & RK(J,I)*U(J,I) GRAD2(J,I)=FN2(J,I)-FM2(J,I)+T2(J,I)/H1(J,I)-FL2(J,I)-YL*U(J,I)- & RK(J,I)*V(J,I) GRAD1(J,I)=GRAD1(J,I)/GRAV*HL1(J,I)*HL1(J,I) GRAD2(J,I)=GRAD2(J,I)/GRAV*HL2(J,I)*HL2(J,I) 1 CONTINUE IF(IBOUND.EQ.1)THEN DO 2 I=1,M DZI(1,I)=0.0 DO 3 J=2,N1 DZI(J,I)=DZI(J-1,I)+GRAD2(J-1,I)*HAG*HL2(J-1,I) 3 CONTINUE 2 CONTINUE ELSE END IF IF(IBOUND.EQ.2)THEN DO 21 I=1,M DZI(N1,I)=0.0 DO 31 J=1,N JJ=N1-J DZI(JJ,I)=DZI(JJ+1,I)-GRAD2(JJ,I)*HAG*HL2(JJ,I) 31 CONTINUE 21 CONTINUE ELSE END IF IF(IBOUND.EQ.3)THEN DO 22 J=1,N DZI(J,1)=0.0 DO 32 I=2,M1 DZI(J,I)=DZI(J,I-1)+GRAD1(J,I-1)*HAG*HL1(J,I-1) 32 CONTINUE 22 CONTINUE ELSE END IF IF(IBOUND.EQ.4)THEN DO 23 J=1,N DZI(J,M1)=0.0 DO 33 I=1,M II=M1-I DZI(J,II)=DZI(J,II+1)-GRAD1(J,II)*HAG*HL1(J,II) 33 CONTINUE 23 CONTINUE ELSE END IF RETURN END SUBROUTINE CHECK c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' HMAX=0. HMIN=1.E10 DO 2 I=1,M1 DO 2 J=1,N1 IF(H(J,I).EQ.PF)GO TO 2 IF(H(J,I).GT.HMAX)HMAX=H(J,I) IF(H(J,I).LT.HMIN)HMIN=H(J,I) 2 CONTINUE IF(HMAX.GT.PF)THEN WRITE(*,*)'Profundidad maxima es mayor que PF' STOP ELSE END IF IF(HMIN.LT.DMIN)THEN WRITE(*,*)'La profundidad minima es menor que DMIN' STOP ELSE END IF DO 12 I=1,M1 DO 12 J=1,N1 HP(J,I)=H(J,I) 12 CONTINUE RETURN END SUBROUTINE ZEROS c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' INCLUDE 'VARS.INC' DO I=1,M1 DO J=1,N1 DELTA(J,I)=0. X0(J,I)=HAG*(I-1)*SCALE+XMIN Y0(J,I)=HAG*(J-1)+YMIN HH1(J,I)=SCALE HH2(J,I)=1. dzi(J,I)=0. BOTANGL(J,I)=0. END DO END DO DO I=1,M DO J=1,N X(J,I)=HAG*(I-0.5)*SCALE+XMIN Y(J,I)=HAG*(J-0.5)+YMIN HL1(J,I)=SCALE HL2(J,I)=1. QSX(J,I)=0. QSY(J,I)=0. QBX(J,I)=0. QBY(J,I)=0. QSXC(J,I)=0. QSYC(J,I)=0. QBXC(J,I)=0. QBYC(J,I)=0. U1NU(J,I)=0. U2NU(J,I)=0. FL1(J,I)=0. FL2(J,I)=0. FM1(J,I)=0. FM2(J,I)=0. FN1(J,I)=0. FN2(J,I)=0. FN1MD(J,I)=0. FN2MD(J,I)=0. RK(J,I)=1.e-04 GRAD1(J,I)=0. GRAD2(J,I)=0. UORB(J,I)=0. UORBF(J,I)=0. UORBB(J,I)=0. FX(J,I)=0. FY(J,I)=0. AFU(J,I)=0. AFV(J,I)=0. U(J,I)=0. V(J,I)=0. UC(J,I)=0. VC(J,I)=0. HW(J,I)=0.000001 FLAMB(J,I)=100000. TAU(J,I)=100000. UCRW(J,I)=0. ! 0 OR A LARGE VALUE ? UR(J,I)=0. UB(J,I)=0. DW(J,I)=0.001 ASYM(J,I)=0. DO KZ=1,10 UV(KZ,J,I)=0. END DO END DO END DO 5 FORMAT(I5,F8.2,2F10.1) IF(IEN.EQ.0)THEN OPEN (1,FILE=LITOP//'CCHANGX.DAT') DO I=1,M1 XCOSTA(I)=0. XUTM1(I)=0. YUTM1(I)=0. WRITE(1,5)I,XCOSTA(I),XUTM1(I),YUTM1(I) END DO CLOSE (1) OPEN (1,FILE=LITOP//'CCHANGY.DAT') DO J=1,N1 YCOSTA(J)=0. XUTM2(J)=0. YUTM2(J)=0. WRITE(1,5)J,YCOSTA(J),XUTM2(J),YUTM2(J) END DO CLOSE (1) ELSE END IF RETURN END SUBROUTINE OUTPUT c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' INCLUDE 'VARS.INC' IF(TYPE1.EQ.2)THEN IF(NPERF1.NE.0)CALL INTERP1 IF(IEN.EQ.0)IEN=1 IEN=IEN+1 OPEN (1,FILE=LITOP//'ENLACE.DAT') WRITE(1,*)IEN CLOSE (1) OPEN (1,FILE=SWANP//'BATINEW.DAT') DO J=1,N0 IF(IFORM.EQ.1)JJ=N0-J+1 ! IF(IFORM.EQ.0)JJ=J ! DO I=1,M0 H2(JJ,I)=H1(JJ,I) IF(H2(JJ,I).EQ.PF)H2(JJ,I)=0. END DO WRITE(1,111)(H2(JJ,I),I=1,M0) END DO 111 FORMAT(1000F10.6) CLOSE (1) OPEN (1,FILE=LITOP//'BOTTOM.DAT') DO J=1,N0+1 IF(IFORM.EQ.1)JJ=(N0+1)-J+1 IF(IFORM.EQ.0)JJ=J write(1,111)(H(JJ,I),I=1,M0+1) END DO CLOSE (1) ELSE END IF CALL CONVERT OPEN (1,FILE=LITOP//'SURF.DAT') OPEN (2,FILE=LITOP//'SURFUV.DAT') OPEN (3,FILE=LITOP//'SURFQS.DAT') OPEN (4,FILE=LITOP//'SURFQB.DAT') OPEN (5,FILE=LITOP//'SEDI.DAT') OPEN (7,FILE=LITOP//'level.DAT') OPEN (8,FILE=LITOP//'level1.DAT') OPEN (9,FILE=LITOP//'CONCEN.DAT') OPEN (10,FILE=LITOP//'UVTOT.DAT') OPEN (11,FILE=LITOP//'GRAFH.DAT') OPEN (16,FILE=LITOP//'waveu.DAT') DO 8 J=1,N1 DO 8 I=1,M1 IF(X0(J,I).EQ.XLIM)GO TO 8 IF(IPRT(1).EQ.1)WRITE(1,*)X0(J,I),Y0(J,I),PSI(J,I) if(h(j,i).eq.pf)go to 8 IF(IPRT(5).EQ.1)WRITE(5,*)X0(J,I),Y0(J,I),DELTA(J,I),H(J,I) 8 continue DO 80 J=1,N WRITE(8,*)J,DZI(J,10) DO 80 I=2,M-1 CCCC IF(H(J,I).EQ.PF.AND.H(J,I+1).EQ.PF)GO TO 80 CCCC IF(GRAD2(J,I).EQ.0.)GO TO 80 IF(X(J,I).EQ.XLIM.OR.Y0(J,I).EQ.YLIM)GO TO 80 IF(IPRT(6).EQ.1)WRITE(7,*)X(J,I),Y0(J,I),DZI(J,I) 80 continue DO 31 I=1,M DO 31 J=1,N IF(X(J,I).EQ.XLIM.OR.Y(J,I).EQ.YLIM)GO TO 31 IF(UR(J,I).EQ.0.)GO TO 809 WRITE(16,*)X(J,I),Y(J,I),UR(J,I),DIR(J,I) 809 CONTINUE UMOD=SQRT(UC(J,I)*UC(J,I)+VC(J,I)*VC(J,I)+1.e-22) ! QS=SQRT(QSXC(J,I)*QSXC(J,I)+QSYC(J,I)*QSYC(J,I)+1.e-22) ! SEE ME QB=SQRT(QBXC(J,I)*QBXC(J,I)+QBYC(J,I)*QBYC(J,I)+1.e-22) ! CALL TANG(UC(J,I),VC(J,I),ANGL) ! CALL TANG(QSXC(J,I),QSYC(J,I),AQS) ! CALL TANG(QBXC(J,I),QBYC(J,I),AQB) ! IF(UMOD.EQ.0.)GO TO 800 IF(IPRT(2).EQ.1)WRITE(2,*)X(J,I),Y(J,I),UMOD,ANGL IF(IPRT(9).EQ.1)WRITE(11,*)X(J,I),Y(J,I),H1(J,I) 800 CONTINUE IF(QS.EQ.0.)GO TO 801 IF(IPRT(3).EQ.1)WRITE(3,*)X(J,I),Y(J,I),QS,AQS 801 IF(QB.EQ.0.)GO TO 802 IF(IPRT(4).EQ.1)WRITE(4,*)X(J,I),Y(J,I),QB,AQB 802 IF(IPRT(7).EQ.1)WRITE(9,*)X(J,I),Y(J,I),CON(2,J,I) 31 CONTINUE DO 32 I=1,M DO 32 J=1,N IF(X(J,I).EQ.XLIM.OR.Y(J,I).EQ.YLIM)GO TO 32 CALL TANG(UC(J,I),VC(J,I),ANGL) IF(UV(1,J,I).EQ.0.)GO TO 32 IF(IPRT(8).EQ.1)WRITE(10,5)X(J,I),Y(J,I),(UV(K,J,I),K=1,10),ANGL 32 CONTINUE 5 FORMAT(13E13.5) CLOSE (1) CLOSE (2) CLOSE (3) CLOSE (4) CLOSE (5) CLOSE (7) CLOSE (8) CLOSE (9) CLOSE (10) CLOSE (11) C------------------------- PROFILE ANALYSIS -------------- OPEN (1,FILE=LITOP//'TALUDX.DAT') OPEN (2,FILE=LITOP//'TALUDY.DAT') IF(NPERF.EQ.0)GO TO 100 NPERFX=0 NPERFY=0 DO 90 II=1,NPERF IF(MX(II).EQ.0)THEN NPERFX=NPERFX+1 DO I=1,M1 XP(NPERFX,I)=X0(NY(II),I) R=H(NY(II),I) IF(R.EQ.PF)R=0. DEEPX(NPERFX,I)=-R R0=H0(NY(II),I) IF(R0.EQ.PF)R0=0. DEEPX0(NPERFX,I)=-R0 END DO C------------------------- DEAN'S PROFILES ----------------- X IF(H0(NY(II),1).EQ.PF)THEN I0=1 IM1=M1 IP=1 ELSE I0=M1 IM1=1 IP=-1 END IF IPE(NPERFX)=IP DIST=0. DO 200 I=I0,IM1,IP IK=I IF(IP.EQ.-1)IK=I0-I+1 IF(H0(NY(II),I).EQ.PF)THEN XP0=XP(NPERFX,I) PDEANX(NPERFX,IK)=0. ELSE DIST=ABS(XP(NPERFX,I)-XP0)+1.E-11 PDEANX(NPERFX,IK)=-DEAN(DIST,ADEAN) END IF 200 CONTINUE C------------------------------------------------------------- ELSE NPERFY=NPERFY+1 DO J=1,N1 YP(J,NPERFY)=Y0(J,MX(II)) R=H(J,MX(II)) IF(R.EQ.PF)R=0. DEEPY(J,NPERFY)=-R R0=H0(J,MX(II)) IF(R0.EQ.PF)R0=0. DEEPY0(J,NPERFY)=-R0 END DO C------------------------- DEAN'S PROFILES ----------------- Y IF(H0(1,MX(II)).EQ.PF)THEN J0=1 JN=N1 JP=1 ELSE J0=N1 JN=1 JP=-1 END IF JPE(NPERFY)=JP DIST=0. DO 201 J=J0,JN,JP JK=J IF(JP.EQ.-1)JK=J0-J+1 IF(H0(J,MX(II)).EQ.PF)THEN YP0=YP(J,NPERFY) PDEANY(JK,NPERFY)=0. ELSE DIST=ABS(YP(J,NPERFY)-YP0)+1.E-11 PDEANY(JK,NPERFY)=-DEAN(DIST,ADEAN) END IF 201 CONTINUE C------------------------------------------------------------- END IF 90 CONTINUE C------------------------------------------------------------- DO II=1,NPERFX IF(IPE(II).EQ.-1)THEN DO I=1,M1 IM1=M1-I+1 PDEANX1(II,I)=PDEANX(II,IM1) END DO DO I=1,M1 PDEANX(II,I)=PDEANX1(II,I) END DO ELSE END IF END DO DO II=1,NPERFY IF(JPE(II).EQ.-1)THEN DO J=1,N1 JN1=N1-J+1 PDEANY1(J,II)=PDEANY(JN1,II) END DO DO J=1,N1 PDEANY(J,II)=PDEANY1(J,II) END DO ELSE END IF END DO C------------------------------------------------------------- C FILMS IF(TYPE1.EQ.1.OR.IPRT(11).EQ.0)GO TO 970 IF(IEN.EQ.2)THEN OPEN (3,FILE=LITOP//'FILMX.DAT',STATUS='NEW') OPEN (4,FILE=LITOP//'FILMY.DAT',STATUS='NEW') TIEMF=1. ELSE OPEN (3,FILE=LITOP//'FILMX.DAT',STATUS='OLD') OPEN (4,FILE=LITOP//'FILMY.DAT',STATUS='OLD') 950 CONTINUE DO I=1,M1 READ(3,11,END=990)TIEMF,(XP1(II,I),PDEANX1(II,I),DEEPX01(II,I), & DEEPX1(II,I),II=1,NPERFX) END DO 990 DO J=1,N1 READ(4,11,END=960)TIEMF,(YP1(J,II),PDEANY1(J,II),DEEPY01(J,II), & DEEPY1(J,II),II=1,NPERFY) END DO GO TO 950 960 TIEMF=TIEMF+1. END IF DO I=1,M1 WRITE(3,11)TIEMF,(XP(II,I),PDEANX(II,I),DEEPX0(II,I), & DEEPX(II,I),II=1,NPERFX) END DO DO J=1,N1 WRITE(4,11)TIEMF,(YP(J,II),PDEANY(J,II),DEEPY0(J,II), & DEEPY(J,II),II=1,NPERFY) END DO CLOSE (3) CLOSE (4) C------------------ END OF FILM 970 CONTINUE IF(IPRT(10).EQ.0)GO TO 100 DO I=1,M1 WRITE(1,11)(XP(II,I),PDEANX(II,I),DEEPX0(II,I), & DEEPX(II,I),II=1,NPERFX) END DO DO J=1,N1 WRITE(2,11)(YP(J,II),PDEANY(J,II),DEEPY0(J,II), & DEEPY(J,II),II=1,NPERFY) END DO 11 FORMAT(1000E17.8) 100 CONTINUE CLOSE (1) CLOSE (2) C-------------------------------------------------- CONTROL FILE FOR MC/NC POINT IF(IPRT(12).EQ.0)GO TO 102 OPEN (1,FILE=LITOP//'CONTROL.DAT') WRITE(1,*)'Coordenadas del punto: ',MC,NC WRITE(1,*)'Profundidad del punto: ',H1(NC,MC) WRITE(1,*)'Velocidad media (m/s): ',U(NC,MC),V(NC,MC), & SQRT(U(NC,MC)*U(NC,MC)+V(NC,MC)*V(NC,MC)+1.e-22) WRITE(1,*)'Corriente de retorno: ',UR(NC,MC) WRITE(1,*)'Corriente de fondo en ola: ',UB(NC,MC) WRITE(1,*)'Altura de ola: ',HW(NC,MC) WRITE(1,*)'Periodo de ola: ',TAU(NC,MC) WRITE(1,*)'Longitud de ola: ',FLAMB(NC,MC) WRITE(1,*)'Velocidad orbital hacia adelante: ',UORBF(NC,MC) WRITE(1,*)'Velocidad orbital hacia atr| s: ',UORBB(NC,MC) WRITE(1,*)'Orbital excursi�n: ',AD(NC,MC) WRITE(1,*)'Espesor de capa l�mite para olas: ',DW(NC,MC) WRITE(1,*)'Direcci�n de ola: ',DIR(NC,MC) CALL TANG(U(NC,MC),V(NC,MC),ANGU) WRITE(1,*)'Direcci�n de corriente: ',ANGU WRITE(1,*)'Diametro de grano D50: ',D50 WRITE(1,*)'Diametro de grano D90: ',D90 WRITE(1,*)'Diametro sedimentol�gico: ',DS WRITE(1,*)'Rugosidad para corrientes: ',Z0 WRITE(1,*)'Rugosidad para oleaje: ',KSW WRITE(1,*)'Factor de fricci�n en olas: ',FW WRITE(1,*)'Factor de eficiencia de corr.: ',FMUC WRITE(1,*)'Factor de eficiencia de olas: ',FMUW WRITE(1,*)'Coeficiente de interacci�n C-O: ',ALFCW WRITE(1,*)'Estr~Bs adim. para transp. de fondo',T(NC,MC) WRITE(1,*)'Estr~Bs adim. para concentracion: ',TA(NC,MC) WRITE(1,*)'Transporte en suspension: ',QSX(NC,MC),QSY(NC,MC), & SQRT(QSX(NC,MC)**2+QSY(NC,MC)**2+1.e-22) CALL TANG(QSX(NC,MC),QSY(NC,MC),ANGU) WRITE(1,*)'Direccion en suspension: ',ANGU WRITE(1,*)'Transporte en el fondo: ',QBX(NC,MC),QBY(NC,MC), & SQRT(QBX(NC,MC)**2+QBY(NC,MC)**2+1.e-22) CALL TANG(QBX(NC,MC),QBY(NC,MC),ANGU) WRITE(1,*)'Direccion en el fondo: ',ANGU CLOSE (1) 102 CONTINUE C------------------------- PUNTOS DE CONTROL H(t) -------- IF(TYPE1.EQ.1.OR.IPRT(13).EQ.0)GO TO 97 IF(IEN.EQ.2)THEN OPEN (1,FILE=LITOP//'HDET.DAT',STATUS='NEW') TIEMP=1. ELSE OPEN (1,FILE=LITOP//'HDET.DAT',STATUS='OLD') 95 READ(1,*,END=96)TIEMP,(HVREM(NCH(K),MCH(K)),K=1,NHC) GO TO 95 96 TIEMP=TIEMP+1. END IF DO K=1,NHC HVREM(NCH(K),MCH(K))=H(NCH(K),MCH(K)) IF(H(NCH(K),MCH(K)).EQ.PF)HVREM(NCH(K),MCH(K))=0. END DO WRITE(1,94)TIEMP,(HVREM(NCH(K),MCH(K)),K=1,NHC) 94 FORMAT(1000F10.3) CLOSE (1) 97 CONTINUE IF(TYPE1.EQ.1)RETURN IF(ICORR.EQ.0)RETURN C------------------------- OLEAJE EN CORRIENTE------------ OPEN(1,FILE=SWANP//'CORIENTE.DAT') DO J=1,N IF(IFORM.EQ.1)JJ=N-J+1 IF(IFORM.EQ.0)JJ=J WRITE(1,101)(U(JJ,I),I=1,M) END DO DO J=1,N IF(IFORM.EQ.1)JJ=N-J+1 IF(IFORM.EQ.0)JJ=J WRITE(1,101)(V(JJ,I),I=1,M) END DO 101 FORMAT(1000F6.2) CLOSE (1) RETURN END REAL FUNCTION DEAN(DIST,ADEAN) DEAN=ADEAN*DIST**(2./3.) RETURN END SUBROUTINE SWITCHE c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' INCLUDE 'VARS.INC' OPEN (1,FILE=LITOP//'PARAM2.DAT') READ(1,*)IWIC READ(1,*)INLT READ(1,*)IWC READ(1,*)IRTF READ(1,*)ITB READ(1,*)ITRANS READ(1,*)ISLOP READ(1,*)DLIMM READ(1,*)DLIMP READ(1,*)QLIM READ(1,*)ICORR READ(1,*)ISMOOT IF(ISMOOT.EQ.1)THEN READ(1,*)KRATS READ(1,*)WEIGHT IF(WEIGHT.GT.1.0.OR.WEIGHT.LT.0.0) & STOP'Error en definicion de WEIGHT' ELSE READ(1,*)KRATS0 READ(1,*)WEIGHT0 END IF READ(1,*)ISMOOTW IF(ISMOOTW.NE.1)GO TO 2 READ(1,*)KRATSW READ(1,*)WEIGHTW IF(WEIGHTW.GT.1.0.OR.WEIGHTW.LT.0.0) & STOP'Error en definicion de WEIGHTW' READ(1,*)ICOMP IF(ICOMP.EQ.2.OR.ICOMP.GT.5)STOP'Error en definicion de ICOMP' 2 CONTINUE IF(IWIC.NE.0.AND.IWIC.NE.1)THEN WRITE(*,*)'Error en definicion de IWIC' stop ELSE END IF IF(INLT.NE.0.AND.INLT.NE.1)THEN WRITE(*,*)'Error en definicion de INLT' stop ELSE END IF IF(IWC.NE.0.AND.IWC.NE.1)THEN WRITE(*,*)'Error en definicion de IWC' stop ELSE END IF IF(IRTF.NE.0.AND.IRTF.NE.1)THEN WRITE(*,*)'Error en definicion de IRTF' stop ELSE END IF IQS=1 IQB=1 IF(ITRANS.EQ.1)IQB=0 IF(ITRANS.EQ.2)IQS=0 CLOSE (1) RETURN END SUBROUTINE CONVERT c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' INCLUDE 'VARS.INC' C---------------- CONVERTIR LOS COMPONENTES CONTRAVARIANTES A LOS CARTESIANOS DO 100 I=2,M-1 DO 100 J=2,N-1 YN=(Y(J+1,I)-Y(J-1,I))*0.5/HAG YK=(Y(J,I+1)-Y(J,I-1))*0.5/HAG XN=(X(J+1,I)-X(J-1,I))*0.5/HAG XK=(X(J,I+1)-X(J,I-1))*0.5/HAG UC(J,I)=U(j,i)*xk+V(j,i)*xn VC(J,I)=U(j,i)*yk+V(j,i)*yn QSXC(J,I)=QSX(j,i)*xk+QSY(j,i)*xn QSYC(J,I)=QSX(j,i)*yk+QSY(j,i)*yn QBXC(J,I)=QBX(j,i)*xk+QBY(j,i)*xn QBYC(J,I)=QBX(j,i)*yk+QBY(j,i)*yn 100 CONTINUE c------------------------------------------------------------------ DO 101 I=2,M-1 J=1 YN=(Y0(J+1,I)-Y0(J,I)+Y0(J+1,I+1)-Y0(J,I+1))*0.5/HAG YK=(Y0(J,I+1)-Y0(J,I)+Y0(J+1,I+1)-Y0(J+1,I))*0.5/HAG XN=(X0(J+1,I)-X0(J,I)+X0(J+1,I+1)-X0(J,I+1))*0.5/HAG XK=(X0(J,I+1)-X0(J,I)+X0(J+1,I+1)-X0(J+1,I))*0.5/HAG UC(J,I)=U(j,i)*xk+V(j,i)*xn VC(J,I)=U(j,i)*yk+V(j,i)*yn QSXC(J,I)=QSX(j,i)*xk+QSY(j,i)*xn QSYC(J,I)=QSX(j,i)*yk+QSY(j,i)*yn QBXC(J,I)=QBX(j,i)*xk+QBY(j,i)*xn QBYC(J,I)=QBX(j,i)*yk+QBY(j,i)*yn 109 J=N YN=(Y0(J+1,I)-Y0(J,I)+Y0(J+1,I+1)-Y0(J,I+1))*0.5/HAG YK=(Y0(J,I+1)-Y0(J,I)+Y0(J+1,I+1)-Y0(J+1,I))*0.5/HAG XN=(X0(J+1,I)-X0(J,I)+X0(J+1,I+1)-X0(J,I+1))*0.5/HAG XK=(X0(J,I+1)-X0(J,I)+X0(J+1,I+1)-X0(J+1,I))*0.5/HAG UC(J,I)=U(j,i)*xk+V(j,i)*xn VC(J,I)=U(j,i)*yk+V(j,i)*yn QSXC(J,I)=QSX(j,i)*xk+QSY(j,i)*xn QSYC(J,I)=QSX(j,i)*yk+QSY(j,i)*yn QBXC(J,I)=QBX(j,i)*xk+QBY(j,i)*xn QBYC(J,I)=QBX(j,i)*yk+QBY(j,i)*yn 101 CONTINUE DO 102 J=2,N-1 I=1 YN=(Y0(J+1,I)-Y0(J,I)+Y0(J+1,I+1)-Y0(J,I+1))*0.5/HAG YK=(Y0(J,I+1)-Y0(J,I)+Y0(J+1,I+1)-Y0(J+1,I))*0.5/HAG XN=(X0(J+1,I)-X0(J,I)+X0(J+1,I+1)-X0(J,I+1))*0.5/HAG XK=(X0(J,I+1)-X0(J,I)+X0(J+1,I+1)-X0(J+1,I))*0.5/HAG UC(J,I)=U(j,i)*xk+V(j,i)*xn VC(J,I)=U(j,i)*yk+V(j,i)*yn QSXC(J,I)=QSX(j,i)*xk+QSY(j,i)*xn QSYC(J,I)=QSX(j,i)*yk+QSY(j,i)*yn QBXC(J,I)=QBX(j,i)*xk+QBY(j,i)*xn QBYC(J,I)=QBX(j,i)*yk+QBY(j,i)*yn 108 I=M YN=(Y0(J+1,I)-Y0(J,I)+Y0(J+1,I+1)-Y0(J,I+1))*0.5/HAG YK=(Y0(J,I+1)-Y0(J,I)+Y0(J+1,I+1)-Y0(J+1,I))*0.5/HAG XN=(X0(J+1,I)-X0(J,I)+X0(J+1,I+1)-X0(J,I+1))*0.5/HAG XK=(X0(J,I+1)-X0(J,I)+X0(J+1,I+1)-X0(J+1,I))*0.5/HAG UC(J,I)=U(j,i)*xk+V(j,i)*xn VC(J,I)=U(j,i)*yk+V(j,i)*yn QSXC(J,I)=QSX(j,i)*xk+QSY(j,i)*xn QSYC(J,I)=QSX(j,i)*yk+QSY(j,i)*yn QBXC(J,I)=QBX(j,i)*xk+QBY(j,i)*xn QBYC(J,I)=QBX(j,i)*yk+QBY(j,i)*yn 102 CONTINUE UC(1,1)=UC(1,2)+UC(2,1)-UC(2,2) UC(N,1)=UC(N,2)+UC(N-1,1)-UC(N-1,2) UC(1,M)=UC(1,M-1)+UC(2,M)-UC(2,M-1) UC(N,M)=UC(N,M-1)+UC(N-1,M)-UC(N-1,M-1) VC(1,1)=VC(1,2)+VC(2,1)-VC(2,2) VC(N,1)=VC(N,2)+VC(N-1,1)-VC(N-1,2) VC(1,M)=VC(1,M-1)+VC(2,M)-VC(2,M-1) VC(N,M)=VC(N,M-1)+VC(N-1,M)-VC(N-1,M-1) QSXC(1,1)=QSXC(1,2)+QSXC(2,1)-QSXC(2,2) QSXC(N,1)=QSXC(N,2)+QSXC(N-1,1)-QSXC(N-1,2) QSXC(1,M)=QSXC(1,M-1)+QSXC(2,M)-QSXC(2,M-1) QSXC(N,M)=QSXC(N,M-1)+QSXC(N-1,M)-QSXC(N-1,M-1) QSYC(1,1)=QSYC(1,2)+QSYC(2,1)-QSYC(2,2) QSYC(N,1)=QSYC(N,2)+QSYC(N-1,1)-QSYC(N-1,2) QSYC(1,M)=QSYC(1,M-1)+QSYC(2,M)-QSYC(2,M-1) QSYC(N,M)=QSYC(N,M-1)+QSYC(N-1,M)-QSYC(N-1,M-1) QBXC(1,1)=QBXC(1,2)+QBXC(2,1)-QBXC(2,2) QBXC(N,1)=QBXC(N,2)+QBXC(N-1,1)-QBXC(N-1,2) QBXC(1,M)=QBXC(1,M-1)+QBXC(2,M)-QBXC(2,M-1) QBXC(N,M)=QBXC(N,M-1)+QBXC(N-1,M)-QBXC(N-1,M-1) QBYC(1,1)=QBYC(1,2)+QBYC(2,1)-QBYC(2,2) QBYC(N,1)=QBYC(N,2)+QBYC(N-1,1)-QBYC(N-1,2) QBYC(1,M)=QBYC(1,M-1)+QBYC(2,M)-QBYC(2,M-1) QBYC(N,M)=QBYC(N,M-1)+QBYC(N-1,M)-QBYC(N-1,M-1) RETURN END SUBROUTINE SHIELD c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 C SHIELD PARAMETER (TETACR) AND CREATICAL BED-SHEAR STRESS (TAUCR) INCLUDE 'COMML.FOR' IF(DCONST.LT.1.)STOP'VERY SMALL PARTICLES OR ANOTHER REASON' IF(DCONST.LE.4.)THEN TETACR=0.24/DCONST TAUCR=(ROS-RO)*GRAV*D50*TETACR RETURN ELSE END IF IF(DCONST.LE.10.)THEN TETACR=0.14/DCONST**0.64 TAUCR=(ROS-RO)*GRAV*D50*TETACR RETURN ELSE END IF IF(DCONST.LE.20.)THEN TETACR=0.04/DCONST**0.1 TAUCR=(ROS-RO)*GRAV*D50*TETACR RETURN ELSE END IF IF(DCONST.LE.150.)THEN TETACR=0.013*DCONST**0.29 TAUCR=(ROS-RO)*GRAV*D50*TETACR RETURN ELSE END IF IF(DCONST.GT.150.)THEN TETACR=0.055 TAUCR=(ROS-RO)*GRAV*D50*TETACR RETURN ELSE END IF RETURN END SUBROUTINE WAVEL c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 C WAVE LENGTH MODIFIED BY CURRENTS & COMPUTED WAVE PARAMETERS INCLUDE 'COMML.FOR' RAD=3.14/180. DO 1 J=1,N DO 1 I=1,M IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 IF(X(J,I).EQ.XLIM)GO TO 1 UMO=SQRT(U(J,I)**2+V(J,I)**2)+0.1E-11 CALL TANG(U(J,I),V(J,I),ANGU) FI1=ABS(ANGU-DIR(J,I)) FLA1(J,I)=FLAMB(J,I) IF(ICORR.EQ.1)GO TO 11 C---------- iterative process ------------------- C KIT=0 C 3 ARG=6.28*H1(J,I)/FLA1(J,I) C KIT=KIT+1 C C IF(UMO.LT.0.01)GO TO 4 C IF(KIT.GT.100000)STOP'ITERATIVE PROCESS DOESN"T CONVERGE FOR LAMB' C TANGH1=(EXP(ARG)-EXP(-ARG))/(EXP(ARG)+EXP(-ARG)) C VREM0=6.28/GRAV*(FLA1(J,I)/TAU(J,I)-UMO*COS(RAD*FI1))**2/TANGH1 C IF(ABS(VREM0-FLA1(J,I)/(VREM0+1.E-11)).LE.0.01)GO TO 2 C IF(VREM0.LT.FLAMB(J,I)/5.)THEN C FLA1(J,I)=FLAMB(J,I)*1.5 C ELSE C FLA1(J,I)=VREM0 C END IF C GO TO 3 C 2 FLAMB(J,I)=VREM0 C------------------------------------------------ ccccc FLAMB(J,I)=FLA1(J,I)-TAU(J,I)*UMO*COS(RAD*FI1) !!!!!!!!!!!!!!! c------------------------------------------------ 4 continue ccccc TAU(J,I)=TAU(J,I)/(1.-UMO*TAU(J,I)*COS(RAD*FI1)/FLAMB(J,I)) !!!! 11 CONTINUE FLAMB(J,I)=AMAX1(FLAMB(J,I),2.0) ! 3.0 m TAU(J,I)=AMAX1(TAU(J,I),0.1) ARG=6.28*H1(J,I)/FLAMB(J,I) AD(J,I)=HW(J,I)/(EXP(ARG)-EXP(-ARG)) if(ad(j,i).eq.0.)write(*,*)i,j,hw(j,i),arg,h1(j,i),flamb(j,i) !!! cccc UORB(J,I)=6.28*HW(J,I)/TAU(J,I)/(EXP(ARG)-EXP(-ARG)) ! why not from SWAN DW(J,I)=0.072*AD(J,I)*(AD(J,I)/KSW)**0.25 C------------------------------------------------- IF(H1(J,I).GE.(0.01*GRAV*TAU(J,I)**2))then ! Tp no modificado ! ? UORBF(J,I)=UORB(J,I)+12.*(6.28*HW(J,I))**2/TAU(J,I)/FLAMB(J,I)/ & (EXP(ARG)-EXP(-ARG))**4 UORBB(J,I)=UORB(J,I)-12.*(6.28*HW(J,I))**2/TAU(J,I)/FLAMB(J,I)/ & (EXP(ARG)-EXP(-ARG))**4 else ALFA1=1.+0.3*(HW(J,I)/H1(J,I)) UORBF(J,I)=ALFA1*UORB(J,I) UORBB(J,I)=(2.-ALFA1)*UORB(J,I) end if C------------------RETURN VELOCITY MASS TRANSPORT & NEAR-BED WAVE-IND. VELOCITY HT=(0.95-0.35*HW(J,I)/H1(J,I))*H1(J,I) UR(J,I)=0.125*SQRT(GRAV/H1(J,I))*HW(J,I)**2/HT ALFAS=UORBF(J,I)/(UORBF(J,I)+UORBB(J,I)+.1e-11) UB(J,I)=(0.05-(ALFAS-0.5))*UORB(J,I) C------------------APPARENT BED ROUGHNESS ----------------- GAMA=0.8+RAD*FI1-0.3*(RAD*FI1)**2 ARGN=GAMA*UORB(J,I)/SQRT(UMO**2+UR(J,I)**2) IF(ARGN.GT.10.)THEN KA(J,I)=10.*Z0 ELSE KA(J,I)=Z0*EXP(ARGN) KA(J,I)=AMIN1(KA(J,I),10.*Z0) IF(KA(J,I).LT.Z0)KA(J,I)=Z0 ! ? artificial END IF C------------------FRICTION FACTORS ----------------------- FCP=0.24/(ALOG10(4.*H1(J,I)/D90))**2 FC=0.24/(ALOG10(12.*H1(J,I)/Z0))**2 FA=0.24/(ALOG10(12.*H1(J,I)/KA(J,I)))**2 FWP=0.3 FW=0.3 IF(AD(J,I).GT.0.001)THEN FWP=EXP(-6.+5.2/(AD(J,I)/3./D90)**0.19) FW=EXP(-6.+5.2/(AD(J,I)/KSW)**0.19) FWP=AMIN1(FWP,0.3) FW=AMIN1(FW,0.3) ELSE END IF C------------------EFFECTIVE TIME-AVERAGED BED-SHEAR STRESS FMUC=FCP/FC FMUW=FWP/FW FMUWA=0.6/DCONST ALFCW=(ALOG(90.*DW(J,I)/KA(J,I))/ALOG(90.*DW(J,I)/Z0))**2* & ((ALOG(30.*H1(J,I)/Z0)-1.)/(ALOG(30.*H1(J,I)/KA(J,I))-1.))**2 ALFCW=AMIN1(ALFCW,1.) TC=RO/8.*FC*SQRT(UMO**2+UR(J,I)**2) TW=RO/4.*FW*UORB(J,I)**2 TCW=TC+TW UPC=SQRT(ALFCW*FMUC*TC/RO) ! Effective bed-shear velocity current T(J,I)=(ALFCW*FMUC*TC+FMUW*TW-TAUCR)/TAUCR T(J,I)=AMAX1(0.,T(J,I)) TA(J,I)=(ALFCW*FMUC*TC+FMUWA*TW-TAUCR)/TAUCR TA(J,I)=AMAX1(0.,TA(J,I)) C----------------- BED-LOAD TRANSPORT ---------------------- ZDELTA=AMAX1(3.*DW(J,I),Z0) VRD=UMO*ALOG(30.*ZDELTA/KA(J,I))/(ALOG(30.*H1(J,I)/KA(J,I))-1.) URD=(UR(J,I)/UMO)*VRD ASYM(J,I)=UORBF(J,I)-UORBB(J,I) SUMUDX=(ASYM(J,I)+UB(J,I)-URD)*COS(RAD*DIR(J,I))+VRD*COS(RAD*ANGU) SUMUDY=(ASYM(J,I)+UB(J,I)-URD)*SIN(RAD*DIR(J,I))+VRD*SIN(RAD*ANGU) SUMUDR=SQRT(SUMUDX**2+SUMUDY**2) ALFA2=ABS(VRD)/(ABS(VRD)+ABS(UORB(J,I))) BETA2=0.25*((ALOG(30.*H1(J,I)/Z0)-1.)/ALOG(30.*ZDELTA/Z0))**2 FCWP=ALFA2*BETA2*FCP+(1.-ALFA2)*FWP TBCWP=0.5*RO*FCWP*SUMUDR**2 GAMA2=1.-SQRT(HW(J,I)/H1(J,I)) GAMA2=AMAX1(0.3,GAMA2) QB=0. CALL SLOPE1(I,J,SUMUDX,SUMUDY,SUMUDR) IF(TBCWP.GT.TAUCR*KB)THEN QB=0.25*GAMA2*ROS*D50/(DCONST**0.3)*SQRT(TBCWP/RO)* & (TBCWP/(TAUCR*KB)-1.)**1.5 QB=QB*AS ELSE END IF QBX(J,I)=(SUMUDX/SUMUDR)*QB QBY(J,I)=(SUMUDY/SUMUDR)*QB 1 CONTINUE RETURN END SUBROUTINE PROFU c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 C VELOCITY DISTRIBUTION OVER TE DEPTH INCLUDE 'COMML.FOR' DO 1 J=1,N DO 1 I=1,M IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 UMO=SQRT(U(J,I)**2+V(J,I)**2) DO 2 K=1,10 Z(K,J,I)=H1(J,I)*(K-1)/9.+AMIN1(Z0,KSW) ! WARNING IF(Z(K,J,I).LT.(3.*DW(J,I)))THEN UDELTA=UMO*ALOG(90.*DW(J,I)/KA(J,I))/ & (ALOG(30.*H1(J,I)/KA(J,I))-1.) UV(K,J,I)=UDELTA*ALOG(30.*Z(K,J,I)/Z0)/ALOG(90.*DW(J,I)/Z0) ELSE UV(K,J,I)=UMO*ALOG(30.*Z(K,J,I)/KA(J,I))/ & (ALOG(30.*H1(J,I)/KA(J,I))-1.) END IF 2 CONTINUE 1 CONTINUE RETURN END SUBROUTINE PROFK c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 C MIXING COEFFICIENT DISTRIBUTION OVER TE DEPTH INCLUDE 'COMML.FOR' DSMIN=0.05 DSMAX=0.2 DO 1 J=1,N DO 1 I=1,M IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 UMO=SQRT(U(J,I)**2+V(J,I)**2) DO 2 K=1,10 EPSC=0. C=18.*ALOG10(12.*H1(J,I)/Z0) UC1=SQRT(GRAV)/C*SQRT(UMO**2+UR(J,I)**2) BETA1=AMIN1(1.5,(1.+2.*(WG/UC1)**2)) IF(Z(K,J,I).LT.(0.5*H1(J,I)))THEN EPSC=0.4*BETA1*UC1*Z(K,J,I)*(1.-Z(K,J,I)/H1(J,I)) ELSE EPSC=0.25*0.4*BETA1*UC1*H1(J,I) END IF EPSW=0. IF(IWIC.EQ.0)GO TO 3 DS0=0.3*H1(J,I)*SQRT(HW(J,I)/H1(J,I)) DS0=AMAX1(DS0,DSMIN) DS0=AMIN1(DS0,DSMAX) EPSBED=0.004*DCONST*DS0*UORB(J,I) EPSTOP=0.035*H1(J,I)*HW(J,I)/TAU(J,I) ! TAU MODIFICADO? IF(Z(K,J,I).LE.DS0)THEN EPSW=EPSBED ELSE IF(Z(K,J,I).GE.(0.5*H1(J,I)))THEN EPSW=EPSTOP ELSE EPSW=EPSBED+(EPSTOP-EPSBED)*(Z(K,J,I)-DS0)/(0.5*H1(J,I)-DS0) END IF 3 CONTINUE EPSCW(K,J,I)=SQRT(EPSC**2+EPSW**2) 2 CONTINUE 1 CONTINUE RETURN END SUBROUTINE FALL c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 C FALL VELOCITY (DS - SIEVE DIAMETER !!!) INCLUDE 'COMML.FOR' IF(DS.LT.1.E-06)THEN WG=0. RETURN ELSE END IF IF(DS.LE.1.E-04)THEN WG=(ROS/RO-1.)*GRAV*DS**2/18./ANU RETURN ELSE END IF IF(DS.LT.1.E-03)THEN WG=10.*ANU/DS*(SQRT(1.+0.01*(ROS/RO-1.)*GRAV*DS**3/ANU**2)-1.) RETURN ELSE WG=1.1*SQRT((ROS/RO-1.)*GRAV*DS) END IF RETURN END SUBROUTINE PROFC c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 C CONCENTRATION DISTRIBUTION OVER TE DEPTH INCLUDE 'COMML.FOR' ALEV=AMAX1(Z0,KSW) ! REFERENCE LEVEL C0=0.65 ! MAXIMUM VOLUME CONCENTRATION DO 1 J=1,N DO 1 I=1,M IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 CBED=0.015*D50/ALEV*TA(J,I)**1.5/DCONST**0.3 CBED=AMIN1(CBED,C0) DO 2 K=1,10 IF(Z(K,J,I).LE.ALEV)THEN CON(K,J,I)=CBED ELSE IF(K.EQ.1)STOP'ERROR IN CONCENTRATION LEVELS' CON(K,J,I)=CON(K-1,J,I)-H1(J,I)/9.* &(1.-CON(K-1,J,I))**5*CON(K-1,J,I)*WG/ &EPSCW(K,J,I)/(1.+(CON(K-1,J,I)/C0)**0.8-2.*(CON(K-1,J,I)/C0)**0.4) CON(K,J,I)=AMAX1(CON(K,J,I),0.0) END IF 2 CONTINUE 1 CONTINUE RETURN END SUBROUTINE FLUX90 c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' RAD=3.14/180. DO 1 J=1,N DO 1 I=1,M QSC=0. QSW=0. IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 DZ=H1(J,I)/9. DO 2 K=1,10 IF(Z(K,J,I).LE.ALEV)GO TO 2 QSC=QSC+ROS*DZ*UV(K,J,I)*CON(K,J,I) QSW=QSW+ROS*DZ*UR(J,I)*CON(K,J,I) 2 CONTINUE CALL TANG(U(J,I),V(J,I),ANGU) QSX(J,I)=QSC*COS(RAD*ANGU)-QSW*COS(RAD*DIR(J,I)) QSY(J,I)=QSC*SIN(RAD*ANGU)-QSW*SIN(RAD*DIR(J,I)) QTX(J,I)=QSX(J,I)*IQS+QBX(J,I)*IQB QTY(J,I)=QSY(J,I)*IQS+QBY(J,I)*IQB 1 CONTINUE RETURN END SUBROUTINE VORTX1 c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' epsp=0.0001 MOL=0 DO 19 I=1,M1 DO 19 J=1,N1 IF(KT.EQ.1)PSI(J,I)=PSI0(J,I) IF(H(J,I).EQ.PF)GOTO 19 MOL=MOL+1 19 CONTINUE TETA=1. L=0 35 L=L+1 IF(L.GE.2000000)GOTO 325 CALL BOUND DO 778 I=2,M DO 778 J=2,N IF(H(J,I).EQ.PF)GOTO 778 C IF(H(J+1,I).EQ.PF.OR.H(J,I+1).EQ.PF.OR.H(J,I-1).EQ.PF C &.OR.H(J-1,I).EQ.PF)GOTO 778 C IF(H(J+1,I+1).EQ.PF.OR.H(J-1,I+1).EQ.PF.OR.H(J+1,I-1).EQ.PF C &.OR.H(J-1,I-1).EQ.PF)GOTO 778 DDT=(AFV(J,I)+AFV(J-1,I)-AFV(J,I-1)-AFV(J-1,I-1) &-AFU(J,I)-AFU(J,I-1)+AFU(J-1,I)+AFU(J-1,I-1)) A1=RK(J,I)*HL1(J,I)/HL2(J,I)/H1(J,I) A4=RK(J,I-1)*HL1(J,I-1)/HL2(J,I-1)/H1(J,I-1) A8=RK(J-1,I)*HL1(J-1,I)/HL2(J-1,I)/H1(J-1,I) A9=RK(J-1,I-1)*HL1(J-1,I-1)/HL2(J-1,I-1)/H1(J-1,I-1) B4=YL*HH1(J+1,I)/HH2(J+1,I)/H(J+1,I) B8=YL*HH1(J-1,I)/HH2(J-1,I)/H(J-1,I) G1=RK(J,I)*HL2(J,I)/HL1(J,I)/H1(J,I) G2=RK(J-1,I)*HL2(J-1,I)/HL1(J-1,I)/H1(J-1,I) G6=RK(J,I-1)*HL2(J,I-1)/HL1(J,I-1)/H1(J,I-1) G7=RK(J-1,I-1)*HL2(J-1,I-1)/HL1(J-1,I-1)/H1(J-1,I-1) O2=YL*HH2(J,I+1)/HH1(J,I+1)/H(J,I+1) O6=YL*HH2(J,I-1)/HH1(J,I-1)/H(J,I-1) PSI1=((G2+G1)*PSI0(J,I+1)+(G7+G6)*PSI(J,I-1)+ & (A4+A1)*PSI0(J+1,I)+(A9+A8)*PSI(J-1,I)+ & 0.5*(B4*(PSI0(J+1,I+1)-PSI0(J+1,I-1))-B8*(PSI0(J-1,I+1)- & PSI(J-1,I-1))-O2*(PSI0(J+1,I+1)-PSI0(J-1,I+1))+ & O6*(PSI0(J+1,I-1)-PSI(J-1,I-1)))-HAG*DDT)/ ! DOESN'T DIVIDE BY 2 AND ONE HAG! & (G1+G2+G6+G7+A1+A4+A8+A9) PSI(J,I)=TETA*PSI1+(1.-TETA)*PSI0(J,I) 778 CONTINUE CALL VELCTY CALL FRIC CALL RHS if(L.eq.1)go to 35 KAP=0 DO 97 I=1,M1 DO 97 J=1,N1 IF(H(J,I).EQ.PF)GOTO 97 IF(PSI0(J,I).EQ.0..AND.PSI(J,I).EQ.0.)GOTO 322 IF(ABS((PSI(J,I)-PSI0(J,I))/(PSI0(J,I)+0.1E-09)).GT.EPSP)GOTO 97 322 KAP=KAP+1 97 CONTINUE DO 2 I=1,M1 DO 2 J=1,N1 IF(H(J,I).EQ.PF)GOTO 2 PSI0(J,I)=PSI(J,I) 2 CONTINUE IF(KAP.EQ.MOL)GOTO 444 IF(L/ITER1*ITER1.EQ.L)WRITE(*,*) &'Acercamiento esta en ',KAP/FLOAT(MOL)*100,' %' IF(KAP/FLOAT(MOL)*100..GE.PERC)GO TO 444 GOTO 35 444 CONTINUE CCC write(*,*)l,kap,mol RETURN 325 WRITE(*,764) 764 FORMAT(2X,'THE NUMBER OF ITERATIONS FOR PSI IS MORE THEN 2e4'/) RETURN END SUBROUTINE TESTWA c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' DO 1 J=1,N DO 1 I=1,M IF(HL(J,I).EQ.PF)GO TO 1 IF(H1(J,I).LT.DMIN)GO TO 1 IF(UORB(J,I).LT.0.)THEN WRITE(*,*)'VELOCIDAD ORBITAL NEGATIVA EN I,J: ',I,J PAUSE ELSE END IF c IF(FX(J,I).EQ.0..AND.FY(J,I).EQ.0.)THEN c WRITE(*,*)'TENSOR DE TENSIONES DE OLAS = 0 EN I,J: ',I,J c PAUSE c ELSE c END IF IF(HW(J,I).LT.0.)THEN WRITE(*,*)'ALTURA DE OLA NEGATIVA EN I,J: ',I,J PAUSE ELSE END IF IF(HW(J,I).LT.0.001)HW(J,I)=0.001 IF(FLAMB(J,I).LT.0.)THEN WRITE(*,*)'LONGITUD DE OLA NEGATIVA EN I,J: ',I,J PAUSE ELSE END IF IF(TAU(J,I).LT.0.)THEN WRITE(*,*)'PERIODO DE OLA NEGATIVO EN I,J: ',I,J PAUSE ELSE END IF 1 CONTINUE DO 10 I=1,M IF(H1(1,I).LT.DMIN)GO TO 10 FX(1,I)=FX(2,I) FY(1,I)=FY(2,I) 10 CONTINUE DO 11 I=1,M IF(H1(N,I).LT.DMIN)GO TO 11 FX(N,I)=FX(N-1,I) FY(N,I)=FY(N-1,I) 11 CONTINUE DO 12 J=1,N IF(H1(J,1).LT.DMIN)GO TO 12 FX(J,1)=FX(J,2) FY(J,1)=FY(J,2) 12 CONTINUE DO 13 J=1,N IF(H1(J,M).LT.DMIN)GO TO 13 FX(J,M)=FX(J,M-1) FY(J,M)=FY(J,M-1) 13 CONTINUE RETURN END SUBROUTINE COSTA c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' C ORIENTACION CUALQUIERA INCLUDE 'VARS.INC' OPEN (1,FILE=LITOP//'CCHANGX.DAT') DO I=1,M1 READ(1,5,END=80)I1,XCOSTA(I),XUTM1(I),YUTM1(I) END DO 80 CLOSE (1) OPEN (1,FILE=LITOP//'CCHANGY.DAT') DO J=1,N1 READ(1,5,END=81)J1,YCOSTA(J),XUTM2(J),YUTM2(J) END DO 81 CLOSE (1) 5 FORMAT(I5,F8.2,2F10.1) OPEN (1,FILE=LITOP//'CCHANGX.DAT') DO 10 I=1,M1 ! NODOS 1-M1 C BUSQUEDA DE LA COSTA DELTAC=0. J=1 IF(H(J,I).EQ.PF)GO TO 32 GO TO 33 32 IF(H(J,I).EQ.PF.AND.H(J+1,I).NE.PF)THEN XUTM1(I)=(X0(J,I)+X0(J-1,I))/2. XUTM1(I)=(X0(J,I)+X0(J-1,I))/2. ! FUE X0(J,I) YUTM1(I)=(Y0(J,I)+Y0(J-1,I))/2. ! IF(ZNAK2(J,I).EQ.PF3)go to 38 DELTAC=DELTA(J+1,I)*HAG*HH2(J+1,I)/H0(J+1,I) GO TO 31 ELSE END IF 38 J=J+1 IF(J.GE.N1)GO TO 31 GO TO 32 33 IF(H(J,I).NE.PF.AND.H(J+1,I).EQ.PF)THEN XUTM1(I)=(X0(J+1,I)+X0(J+2,I))/2. ! FUE X0(J+1,I) YUTM1(I)=(Y0(J+1,I)+Y0(J+2,I))/2. ! IF(ZNAK2(J+1,I).EQ.PF3)go to 39 DELTAC=DELTA(J,I)*HAG*HH2(J,I)/H0(J,I) GO TO 31 ELSE END IF 39 J=J+1 IF(J.GE.N1)GO TO 10 GO TO 33 31 CONTINUE XCOSTA(I)=XCOSTA(I)+DELTAC IF(IPRT(14).EQ.1)WRITE(1,5)I,XCOSTA(I),XUTM1(I),YUTM1(I) 10 CONTINUE CLOSE (1) C---------------------- COSTA MERIDIANAL -------------- OPEN (1,FILE=LITOP//'CCHANGY.DAT') DO 1 J=1,N1 ! NODOS 1-M1 C BUSQUEDA DE LA COSTA DELTAC=0. I=1 IF(H(J,I).EQ.PF)GO TO 320 GO TO 330 320 IF(H(J,I).EQ.PF.AND.H(J,I+1).NE.PF)THEN XUTM2(J)=(X0(J,I)+X0(J,I-1))/2. ! FUE X0(J,I) YUTM2(J)=(Y0(J,I)+Y0(J,I-1))/2. IF(ZNAK2(J,I).EQ.PF3)go to 380 DELTAC=DELTA(J,I+1)*HAG*HH1(J,I+1)/H0(J,I+1) GO TO 310 ELSE END IF 380 I=I+1 IF(I.GE.M1)GO TO 310 GO TO 320 330 IF(H(J,I).NE.PF.AND.H(J,I+1).EQ.PF)THEN XUTM2(J)=(X0(J,I+1)+X0(J,I+2))/2. ! FUE X0(J,I+1) YUTM2(J)=(Y0(J,I+1)+Y0(J,I+2))/2. ! IF(ZNAK2(J,I+1).EQ.PF3)go to 390 DELTAC=DELTA(J,I)*HAG*HH1(J,I)/H0(J,I) GO TO 310 ELSE END IF 390 I=I+1 IF(I.GE.M1)GO TO 1 GO TO 330 310 CONTINUE YCOSTA(J)=YCOSTA(J)+DELTAC IF(IPRT(14).EQ.1)WRITE(1,5)J,YCOSTA(J),XUTM2(J),YUTM2(J) 1 CONTINUE CLOSE (1) RETURN END SUBROUTINE LIMIT_S c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' DO 1 I=1,M DO 1 J=1,N QMOD=SQRT(QTX(J,I)**2+QTY(J,I)**2) IF(QMOD.EQ.0.)GO TO 1 QMOD1=AMIN1(QMOD,QLIM) QTX(J,I)=QTX(J,I)*QMOD1/QMOD QTY(J,I)=QTY(J,I)*QMOD1/QMOD 1 CONTINUE RETURN END SUBROUTINE SLOPE c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2000 INCLUDE 'COMML.FOR' INCLUDE 'VARS.INC' DO 1 I=2,M DO 1 I=2,M DO 1 J=2,N IF(H(J,I).EQ.PF)GO TO 1 HOP1=H(J,I+1) IF(HOP1.EQ.PF)HOP1=0. HOP2=H(J,I-1) IF(HOP2.EQ.PF)HOP2=0. HOP3=H(J+1,I) IF(HOP3.EQ.PF)HOP3=0. HOP4=H(J-1,I) IF(HOP4.EQ.PF)HOP4=0. ANGX=0. ANGY=0. IF(X0(J,I+1).EQ.XLIM)GO TO 2 IF(X0(J,I-1).EQ.XLIM)GO TO 2 IF(IXYP.EQ.1)ANGX=(HOP1-HOP2)/(X0(J,I+1)-X0(J,I-1)) IF(IXYP.EQ.2)ANGX=(HOP1-HOP2)/(Y0(J,I+1)-Y0(J,I-1)) 2 IF(Y0(J+1,I).EQ.YLIM)GO TO 3 IF(Y0(J-1,I).EQ.YLIM)GO TO 3 IF(IXYP.EQ.1)ANGY=(HOP3-HOP4)/(Y0(J+1,I)-Y0(J-1,I)) IF(IXYP.EQ.2)ANGY=(HOP3-HOP4)/(X0(J+1,I)-X0(J-1,I)) 3 BOTANGL(J,I)=AMAX1(abs(ANGX),abs(ANGY)) 1 CONTINUE OPEN (13,FILE=LITOP//'SLOPE.DAT') DO 4 J=1,N1 DO 4 I=1,M1 IF(X0(J,I).EQ.XLIM.OR.Y0(J,I).EQ.YLIM)GO TO 4 WRITE(13,*)X0(J,I),Y0(J,I),BOTANGL(J,I) 4 CONTINUE CLOSE (13) RETURN END SUBROUTINE SLOPE1(I,J,SUMUDX,SUMUDY,SUMUDR) c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2001 INCLUDE 'COMML.FOR' KB=1. ANGX=0. ANGY=0. IF(ISLOP.EQ.0)RETURN IP=MIN0(I+1,M) IM1=MAX0(I-1,1) JP=MIN0(J+1,N) JM1=MAX0(J-1,1) HOP1=H1(J,IP) IF(HL(J,IP).EQ.PF)HOP1=0. HOP2=H1(J,IM1) IF(HL(J,IM1).EQ.PF)HOP2=0. HOP3=H1(JP,I) IF(HL(JP,I).EQ.PF)HOP3=0. HOP4=H1(JM1,I) IF(HL(JM1,I).EQ.PF)HOP4=0. IF(X(J,IP).EQ.XLIM)GO TO 2 IF(X(J,IM1).EQ.XLIM)GO TO 2 IF(IXYP.EQ.1)ANGX=(HOP1-HOP2)/(X(J,IP)-X(J,IM1)) IF(IXYP.EQ.2)ANGX=(HOP1-HOP2)/(Y(J,IP)-Y(J,IM1)) 2 IF(Y(JP,I).EQ.YLIM)GO TO 3 IF(Y(JM1,I).EQ.YLIM)GO TO 3 IF(IXYP.EQ.1)ANGY=(HOP3-HOP4)/(Y(JP,I)-Y(JM1,I)) IF(IXYP.EQ.2)ANGY=(HOP3-HOP4)/(X(JP,I)-X(JM1,I)) 3 CONTINUE BETAF=-(SUMUDX*ANGX+SUMUDY*ANGY)/SUMUDR KB=SIN(TANTETA+BETAF)/SIN(TANTETA) AS=ATAN(TANTETA)/COS(BETAF)/(ATAN(TANTETA)+ATAN(BETAF))*2. ! 2. -artific. parameter (adjustment) RETURN END SUBROUTINE BOUND c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.1, 2001 INCLUDE 'COMML.FOR' c---------------- eastern and western bounds ----- DO 1 J=2,N if(h(j,M1).eq.pf)go to 1 IF(ITB.EQ.1)THEN PSI(J,M1)=PSI(j,m)*2.-psi(j,m-1) PSI0(J,M1)=PSI(j,m)*2.-psi(j,m-1) ELSE PSI(J,M1)=PSI(j,m) PSI0(J,M1)=PSI(j,m) END IF 1 continue DO 4 J=2,N if(h(j,1).eq.pf)go to 4 IF(ITB.EQ.1)THEN PSI(J,1)=PSI(j,2)*2.-psi(j,3) PSI0(J,1)=PSI(j,2)*2.-psi(j,3) ELSE PSI(J,1)=PSI(j,2) PSI0(J,1)=PSI(j,2) END IF 4 continue c---------------- southern-nothern bounds ----- DO 3 i=2,M if(h(1,i).eq.pf)go to 2 IF(ITB.EQ.1)THEN PSI(1,I)=PSI(2,I)*2.-PSI(3,I) PSI0(1,I)=PSI(2,I)*2.-PSI(3,I) ELSE PSI(1,I)=PSI(2,I) PSI0(1,I)=PSI(2,I) END IF 2 continue if(h(N1,i).eq.pf)go to 3 IF(ITB.EQ.1)THEN PSI(N1,I)=PSI(N,I)*2.-PSI(N-1,I) PSI0(N1,I)=PSI(N,I)*2.-PSI(N-1,I) ELSE PSI(N1,I)=PSI(N,I) PSI0(N1,I)=PSI(N,I) END IF 3 continue C--------------- COSMETICS -------------------- IF(H(1,M1).NE.PF)PSI(1,M1)=PSI(1,M)+PSI(2,M1)-PSI(2,M) IF(H(N1,M1).NE.PF)PSI(N1,M1)=PSI(N1,M)+PSI(N,M1)-PSI(N,M) IF(H(N1,1).NE.PF)PSI(N1,1)=PSI(N1,2)+PSI(N,1)-PSI(N,2) IF(H(1,1).NE.PF)PSI(1,1)=PSI(1,2)+PSI(2,1)-PSI(2,2) IF(H(1,M1).NE.PF)PSI0(1,M1)=PSI(1,M)+PSI(2,M1)-PSI(2,M) IF(H(N1,M1).NE.PF)PSI0(N1,M1)=PSI(N1,M)+PSI(N,M1)-PSI(N,M) IF(H(N1,1).NE.PF)PSI0(N1,1)=PSI(N1,2)+PSI(N,1)-PSI(N,2) IF(H(1,1).NE.PF)PSI0(1,1)=PSI(1,2)+PSI(2,1)-PSI(2,2) C------------------------------------------------- IF(UNI.EQ.0.)RETURN DO 10 I=2,M1 DO 10 J=2,N1 IF(ZNAK1(J,I).EQ.PF1.AND.ZNAK1(J-1,I).NE.PF1)VALOR=PSI(J-1,I) 10 CONTINUE DO 11 I=2,M1 DO 11 J=2,N1 IF(ZNAK1(J,I).EQ.PF1)PSI(J,I)=VALOR 11 CONTINUE RETURN END SUBROUTINE SMOOTH c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' REAL HS(400,400),HS1(400,400) IF(M1.GT.400)STOP'Cambiar M1 en SMOOTH' IF(N1.GT.400)STOP'Cambiar N1 en SMOOTH' DO I=1,M1 DO J=1,N1 HS(J,I)=H(J,I) IF(HS(J,I).EQ.PF)HS(J,I)=0. END DO END DO DO 4 I=2,M DO 4 J=2,N IF(HS(J,I).eq.0.)go to 4 HS1(J,I)=HS(J,I) IF(ZNAK2(J+1,I).EQ.PF3.OR.ZNAK2(J-1,I).EQ.PF3.OR.ZNAK2(J,I-1). & EQ.PF3.OR. & ZNAK2(J,I+1).EQ.PF3.OR.ZNAK2(J+1,I+1).EQ.PF3.OR.ZNAK2(J+1,I-1). & EQ.PF3.OR. & ZNAK2(J-1,I+1).EQ.PF3.OR.ZNAK2(J-1,I-1).EQ.PF3)GO TO 4 IF(H(J+1,I).EQ.PF.OR.H(J-1,I).EQ.PF.OR.H(J,I-1). ! & EQ.PF.OR. ! & H(J,I+1).EQ.PF.OR.H(J+1,I+1).EQ.PF.OR.H(J+1,I-1). ! VSTAVKA !!! & EQ.PF.OR. ! & H(J-1,I+1).EQ.PF.OR.H(J-1,I-1).EQ.PF)GO TO 4 ! HS1(J,I)=0.4*HS(J,I)+0.1*(HS(J+1,I)+HS(J-1,I)+HS(J,I-1)+ & HS(J,I+1))+0.05*(HS(J+1,I+1)+HS(J-1,I+1)+HS(J+1,I-1)+ & HS(J-1,I-1)) HS1(J,I)=AMAX1(HS1(J,I),DMIN) 4 CONTINUE DO 5 I=2,M DO 5 J=2,N IF(H(J,I).EQ.PF)GO TO 5 H(J,I)=HS1(J,I)*WEIGHT+(1.-WEIGHT)*H(J,I) 5 CONTINUE RETURN END subroutine inic include 'VARS.INC' c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 open (1,file=LITOP//'monitor.dat') write(1,*)'0' close (1) return end SUBROUTINE INTERP1 c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' INCLUDE 'VARS.INC' REAL X00(160000),Y00(160000),HI(160000),R(160000),XN(160000) & ,YN(160000),HN(160000),A(160000) IF(IAUTO.EQ.0)THEN RMIN0=HAG RMIN1=HAG L5=5 ELSE END IF L0=0 DO I=1,M DO J=1,N L0=L0+1 X00(L0)=X(J,I) Y00(L0)=Y(J,I) HI(L0)=H1(J,I) END DO END DO IF(L0.GT.160000)STOP'Cambiar dimensiones en INTERP1' C------------------------------------------------------ DO 33 II=1,NPERF1 DO 3 JJ=1,NP(II) rmin=rmin0 XI=XPP(II,JJ) YI=YPP(II,JJ) c ---- Determination of the weights -------- 800 LL=0 DO 20 k=1,L0 R0=SQRT((XI-X00(k))*(XI-X00(k))+(YI-Y00(k))*(YI-Y00(k))) IF(R0.GT.RMIN)GO TO 20 LL=LL+1 R(LL)=R0 XN(LL)=X00(K) YN(LL)=Y00(K) HN(LL)=HI(K) 20 continue IF(LL.LT.L5)THEN RMIN=RMIN+RMIN1 GO TO 800 ELSE END IF NNN=LL j=0 do 22 k1=1,NNN do 22 k=1,NNN j=j+1 A(j)=Sqrt((Xn(k1)-Xn(k))*(Xn(k1)-Xn(k))+(Yn(k1)-Yn(k))* *(Yn(k1)-Yn(k))) 22 continue CALL GELG(R,A,NNN,1,.1e-06,IER) R1=0. R2=0. DO 2 k=1,NNN R1=R1+R(k)*HN(k) R2=R2+R(k) 2 continue HPERF(II,JJ)=-R1/R2 3 continue 33 continue c---------------------------------------------------------- OPEN (1,FILE=LITOP//'PERFILD.DAT') WRITE(1,*)NPERF1 DO II=1,NPERF1 WRITE(1,*)NP(II) DO JJ=1,NP(II) WRITE(1,*)XPP(II,JJ),YPP(II,JJ),HPERF(II,JJ) END DO END DO CLOSE (1) c---------------------------------------------------------- if(IER.eq.-1)write(*,*) *'NO RESULT BECAUSE OF M LESS THAN 1 OR PIVOT ELEMENT AT ANY', *'ELIMINATION STEP EQUAL TO 0' if(ier.eq.0.or.ier.eq.-1)go to 110 write(*,*)'IER=K= ',IER,'WARNING DUE TO POSSIBLE LOSS OF', *'SIGNIFICANCE INDICATED AT ELIMINATION STEP K+1, WHERE PILOT', *'ELEMENT WAS LESS THAN OR EQUAL TO THE INTERNAL TOLERANCE EPS', *'TIMES ABSOLUTELY GREATEST ELEMENT OF MATRIX A.' PAUSE 110 continue RETURN END C .................................................................. C C SUBROUTINE GELG C C PURPOSE C TO SOLVE A GENERAL SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS. C C USAGE C CALL GELG(R,A,M,N,EPS,IER) C C DESCRIPTION OF PARAMETERS C R - THE M BY N MATRIX OF RIGHT HAND SIDES. (DESTROYED) C ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS. C A - THE M BY M COEFFICIENT MATRIX. (DESTROYED) C M - THE NUMBER OF EQUATIONS IN THE SYSTEM. C N - THE NUMBER OF RIGHT HAND SIDE VECTORS. C EPS - AN INPUT CONSTANT WHICH IS USED AS RELATIVE C TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE. C IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS C IER=0 - NO ERROR, C IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR C PIVOT ELEMENT AT ANY ELIMINATION STEP C EQUAL TO 0, C IER=K - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI- C CANCE INDICATED AT ELIMINATION STEP K+1, C WHERE PIVOT ELEMENT WAS LESS THAN OR C EQUAL TO THE INTERNAL TOLERANCE EPS TIMES C ABSOLUTELY GREATEST ELEMENT OF MATRIX A. C C REMARKS C INPUT MATRICES R AND A ARE ASSUMED TO BE STORED COLUMNWISE C IN M*N RESP. M*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN C SOLUTION MATRIX R IS STORED COLUMNWISE TOO. C THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS C GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS C ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN - C INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL C SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE C INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS C GIVEN IN CASE M=1. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH C COMPLETE PIVOTING. C C .................................................................. C SUBROUTINE GELG(R,A,M,N,EPS,IER) INCLUDE 'VARS.INC' C C dimension R(1),A(1) IF(M)23,23,1 C C SEARCH FOR GREATEST ELEMENT IN MATRIX A 1 IER=0 PIV=0. MM=M*M NM=N*M DO 3 L=1,MM TB=ABS(A(L)) IF(TB-PIV)3,3,2 2 PIV=TB I=L 3 CONTINUE TOL=EPS*PIV C A(I) IS PIVOT ELEMENT. PIV CONTAINS THE ABSOLUTE VALUE OF A(I). C C C START ELIMINATION LOOP LST=1 DO 17 K=1,M C C TEST ON SINGULARITY IF(PIV)23,23,4 4 IF(IER)7,5,7 5 IF(PIV-TOL)6,6,7 6 IER=K-1 7 PIVI=1./A(I) J=(I-1)/M I=I-J*M-K J=J+1-K C I+K IS ROW-INDEX, J+K COLUMN-INDEX OF PIVOT ELEMENT C C PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R DO 8 L=K,NM,M LL=L+I TB=PIVI*R(LL) R(LL)=R(L) 8 R(L)=TB C C IS ELIMINATION TERMINATED IF(K-M)9,18,18 C C COLUMN INTERCHANGE IN MATRIX A 9 LEND=LST+M-K IF(J)12,12,10 10 II=J*M DO 11 L=LST,LEND TB=A(L) LL=L+II A(L)=A(LL) 11 A(LL)=TB C C ROW INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A 12 DO 13 L=LST,MM,M LL=L+I TB=PIVI*A(LL) A(LL)=A(L) 13 A(L)=TB C C SAVE COLUMN INTERCHANGE INFORMATION A(LST)=J C C ELEMENT REDUCTION AND NEXT PIVOT SEARCH PIV=0. LST=LST+1 J=0 DO 16 II=LST,LEND PIVI=-A(II) IST=II+M J=J+1 DO 15 L=IST,MM,M LL=L-J A(L)=A(L)+PIVI*A(LL) TB=ABS(A(L)) IF(TB-PIV)15,15,14 14 PIV=TB I=L 15 CONTINUE DO 16 L=K,NM,M LL=L+J 16 R(LL)=R(LL)+PIVI*R(L) 17 LST=LST+M C END OF ELIMINATION LOOP C BACK SUBSTITUTION AND BACK INTERCHANGE 18 IF(M-1)23,22,19 19 IST=MM+M LST=M+1 DO 21 I=2,M II=LST-I IST=IST-LST L=IST-M L=A(L)+.5 DO 21 J=II,NM,M TB=R(J) LL=J DO 20 K=IST,MM,M LL=LL+1 20 TB=TB-A(K)*R(LL) K=J+L R(J)=R(K) 21 R(K)=TB 22 RETURN C ERROR RETURN 23 IER=-1 RETURN END SUBROUTINE SMOOTHW c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' REAL HS(400,400),HS1(400,400) IF(M.GT.399)STOP'Cambiar M en SMOOTHW' IF(N.GT.399)STOP'Cambiar N en SMOOTHW' WEIGH1=WEIGHTW ! PESOS DO 1 LTIME=1,KRATSW ! VECES DO K=1,ICOMP ! COMPONENTES DO I=1,M DO J=1,N IF(K.EQ.1)HS(J,I)=UORB(J,I) IF(K.EQ.2)HS(J,I)=FX(J,I) IF(K.EQ.3)HS(J,I)=FY(J,I) IF(K.EQ.4)HS(J,I)=HW(J,I) IF(K.EQ.5)HS(J,I)=FLAMB(J,I) IF(H1(J,I).EQ.0.)HS(J,I)=0. END DO END DO DO 4 I=2,M-1 DO 4 J=2,N-1 IF(H1(J,I).eq.0.)go to 4 HS1(J,I)=HS(J,I) C IF(H1(J+1,I).EQ.0..OR.H1(J-1,I).EQ.0..OR.H1(J,I-1). ! C & EQ.0..OR. ! C & H1(J,I+1).EQ.0..OR.H1(J+1,I+1).EQ.0..OR.H1(J+1,I-1). ! VSTAVKA !!! C & EQ.0..OR. ! C & H1(J-1,I+1).EQ.0..OR.H1(J-1,I-1).EQ.0.)GO TO 4 ! HS1(J,I)=0.4*HS(J,I)+0.1*(HS(J+1,I)+HS(J-1,I)+HS(J,I-1)+ & HS(J,I+1))+0.05*(HS(J+1,I+1)+HS(J-1,I+1)+HS(J+1,I-1)+ & HS(J-1,I-1)) 4 CONTINUE DO I=2,M IF(H1(1,I).NE.0.)HS1(1,I)=HS1(2,I) IF(H1(N,I).NE.0.)HS1(N,I)=HS1(N-1,I) END DO DO J=2,N IF(H1(J,1).NE.0.)HS1(J,1)=HS1(J,2) IF(H1(J,M).NE.0.)HS1(J,M)=HS1(J,M-1) END DO IF(H1(1,1).NE.0.)HS1(1,1)=HS1(1,2)+HS1(2,1)-HS1(2,2) IF(H1(1,M).NE.0.)HS1(1,M)=HS1(1,M-1)+HS1(2,M)-HS1(2,M-1) IF(H1(N,1).NE.0.)HS1(N,1)=HS1(N-1,1)+HS1(N,2)-HS1(N-1,2) IF(H1(N,M).NE.0.)HS1(N,M)=HS1(N-1,M)+HS1(N,M-1)-HS1(N-1,M-1) DO 5 I=2,M DO 5 J=2,N IF(H1(J,I).EQ.0.)GO TO 5 IF(K.EQ.1)UORB(J,I)=HS1(J,I)*WEIGH1+UORB(J,I)*(1.-WEIGH1) IF(K.EQ.2)FX(J,I)=HS1(J,I)*WEIGH1+FX(J,I)*(1.-WEIGH1) IF(K.EQ.3)FY(J,I)=HS1(J,I)*WEIGH1+FY(J,I)*(1.-WEIGH1) IF(K.EQ.4)HW(J,I)=HS1(J,I)*WEIGH1+HW(J,I)*(1.-WEIGH1) IF(K.EQ.5)FLAMB(J,I)=HS1(J,I)*WEIGH1+FLAMB(J,I)*(1.-WEIGH1) 5 CONTINUE END DO 1 CONTINUE RETURN END SUBROUTINE PRNH c Elaborado por Serguei A. Lonin, 1999-2000 Version 2.2, 2001 INCLUDE 'COMML.FOR' INCLUDE 'VARS.INC' OPEN (1,FILE=LITOP//'BOT0.DAT') DO 8 J=1,N1 DO 8 I=1,M1 IF(X0(J,I).EQ.XLIM)GO TO 8 if(h0(j,i).eq.pf)go to 8 IF(IPRT(5).EQ.1)WRITE(1,*)X0(J,I),Y0(J,I),H0(J,I) 8 continue RETURN END