lito.for

     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