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
Comments
0 B
|👍
/👎
0 B
|👍
/👎