From 93afe9a7d61dd544b32c5abf3f6dc8a46e57a943 Mon Sep 17 00:00:00 2001
From: Theophile Terraz <theophile.terraz@inrae.fr>
Date: Thu, 9 Dec 2021 15:09:03 +0100
Subject: [PATCH] sairube V091116
---
sairube.f | 5180 +++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 5180 insertions(+)
create mode 100644 sairube.f
diff --git a/sairube.f b/sairube.f
new file mode 100644
index 0000000..f92fe3d
--- /dev/null
+++ b/sairube.f
@@ -0,0 +1,5180 @@
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+C C
+C PROGRAMME DE SAISIE DES DONNEES ASSOCIEES AU PROGRAMME C
+C DE PROPAGATION D'ONDE DE RUPTURE DE BARRAGE C
+C C
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+C C
+C VERSION INITIALE DU 26/09/1988 AUTEUR : JP VILA C
+C MODIFICATIONS PAR G.MARTINET ET A.PAQUIER C
+C VERSION DU 19/6/98 DE A.PAQUIER C
+C version du 15 octobre 2007 :ecriture geometrie dans geomac-i
+C corrige le 22/10
+C corrections du 17/11/08 sur :
+c les parametres des ouvrages y compris largeur maximale de B
+C et le controle de l'ordre des abscisses
+C mise a jour octobre 2009 pour rubarBE
+C le 25 aout 2016 ajout des points 10 et 11 du menu
+C pour passer de monbief a multibief
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+C C
+C LISTE FICHIERS C
+C ---------------- C
+C C
+C DONNEE.ETUDE : CONTIENT PARAMETRES DE CALCUL C
+C ts.etude : PARAMETRES DU TRANSPORT SOLIDE C
+C MAIL0 : ABSCISSES MAILLAGE ET PAS D ESPACE C
+C GEOMAC-I00 : GEOMETRIE C
+C geomac-i0 : TOPOGRAPHIE ABSCISSE-COTE (AUX INTERMAILLES) C
+C geomac-i : GEOMETRIE (TOPO+SÉD) ABSCISSE-COTE (AUX INTERMAILLES) C
+C frot0 : DONNEES FROTTEMENTS C
+C frot : FROTTEMENTS (AUX CENTREMAILLES) C
+C condi0 : DONNEES CONDITIONS INITIALES HYDRAULIQUES C
+C condin : CONDITIONS INITIALES HYDRAULIQUES (AUX CENTREMAILLES) C
+C condis0 : DONNEES CONDITIONS INITIALES SEDIMENTAIRES C
+C condins : CONDITIONS INITIALES SEDIMENTAIRES (AUX CENTREMAILLES)C
+C ABSHYD : ABSCISSE DES HYDROGRAMMES C
+C TNPROF : TEMPS DE SAUVEGARDE DES PROFILS C
+C C
+C LES FICHIERS'0' INDIQUENT QU'ILS CONTIENNENT LES C
+C VARIABLES INTRODUITES AU CLAVIER ET CEUX SANS LE C
+C '0' QU'ILS CONTIENNENT LES VARIABLES INTERPOLEES C
+C C
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+C C
+C ETUDE EST LE NOM DU TRAVAIL.L'UTILISATION DU CODE C
+C SE DEROULE EN PLUSIEURS ETAPES C
+C C
+C 1- GENERATION DU MAILLAGE FICHIER MAIL.ETUDE C
+C 2- CREATION DES FICHIERS GEOM0.ETUDE,FROT0.ETUDE C
+C ,CONDI0.ETUDE EN SELECTIONNANT UN CERTAIN C
+C NOMBRE DE SECTIONS CARACTERISTIQUES LE LONG C
+C DU PROFIL C
+C 3- CREATION DES FICHIERS GEOMac-i.ETUDE,FROT.ETUDE C
+C ,CONDIN.ETUDE PAR INTERPOLATION LINEAIRE SUR C
+C LE MAILLAGE C
+C---------------------------------------------------------C
+C STRUCTURE DES ET DENOMINATION DES FICHIERS: C
+C C
+C LM EST LE NOMBRE DE POINTS DE CALCULS C
+C C
+C-->MAILLAGE : C
+C------------- C
+C MAIL0.ETUDE C
+C BOUCLE SUR LES SECTIONS DE DONNEE C
+C X, DX (PAS D ESPACE DESIRE) C
+C FIN DE LA BOUCLE C
+C MAIL.ETUDE C
+C LM = NOMBRE DE POINTS DE CALCUL C
+C BOUCLE SUR I DE 1 A LM SUR ABSCISSES DES POINTS C
+C DE CALCUL : TMAIL(I) PUIS SUR XTMAIL(I) C
+C ABSCISSES DES INTERMAILLES C
+C-->GEOMETRIE : C
+C-------------- C
+C geomac-i0.etude C
+C LM-1 [nb de sections du fichier] C
+C Pour chaque section: C
+C I,XSECT(I),IJ0(I) [num section, abscisse section, nb couples]C
+C Pour chaque couple de points de la section: C
+C CCOU(J),YCOU(J),ZCOU(J) [marque, abscisse latérale, cote] C
+C C
+C geomac-i.etude C
+C LM-1 [nb de sections du fichier] C
+C Pour chaque section: C
+C I,XSECT(I),IJ0(I) [num section, abscisse section, nb couples]C
+C Pour chaque couple de points de la section: C
+C CCOU(J),YCOU(J),+NBCS(J)fois(ZCS(J,K),DCS(J,K),SCS(J,K),TMCS(J,K))
+C [marque, abscisse latérale, puis pour chaque compart séd:] C
+C [(cote du toit, diamètre, étendue, contrainte mise en mvt)]C
+C C
+C-->FROTTEMENTS : C
+C---------------- C
+C FROT0.ETUDE C
+C BOUCLE SUR LES SECTIONS DE DONNEE C
+C X, FROT0 (COEFFICIENT DE STRICKLER) C
+C C
+C FROT.ETUDE C
+C BOUCLE SUR I DE 1 A LM C
+C I, FROT(I) (COEFFICIENT DE STRICKLER) C
+C C
+C-->CONDITIONS INITIALES : C
+C------------------------- C
+C CONDI0.ETUDE C
+C BOUCLE SUR LES SECTIONS DE DONNEE C
+C X, Y (TIRANT D'EAU), V (VITESSE),H(COTE) C
+C C
+C CONDIN.ETUDE C
+C BOUCLE SUR I DE 1 A LM C
+C I, Y(I) (TIRANT D'EAU), V(I) (VITESSE) C
+C C
+C-->CONDITIONS AUX LIMITES : C
+C--------------------------- C
+C CONDAV.ETUDE C
+C QYH(I), YH(I) C
+C C
+C HYDRO.ETUDE C
+C TH(I), QTH(I) C
+C C
+C HAMONT.ETUDE C
+C HIMP(I),TA(I) C
+C C
+C-->AUTRES CAS : C
+C--------------- C
+C PARAMETRES DE SAUVEGARDES : TEMPS ET ABSCISSES C
+C " " DEVERSEMENT C
+C VISUALISATION DES VARIABLES C
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
+ PROGRAM SAISIE
+
+C IMPLICIT NONE
+
+ INTEGER LMAX,NCMAX,NSMAX,IDPARS,IE,NCOURB
+ PARAMETER (LMAX=3000,NCMAX=300,NSMAX=3000)
+ PARAMETER(NCOURB=10)
+
+ DOUBLE PRECISION eps
+C :CFL,
+c + TINIT,
+c + TMAX,
+c + TIOPDT,
+c + DT,
+c + TS,
+c + PSAVE,DTSAUV,
+c + FDEB1,FDEB2,FDEB3,EPS
+
+ CHARACTER ETUDE*20
+ + ,NOMFIC*38
+ + ,IODEV*1
+ + ,IOPDT*1
+ + ,IOVIS*1
+ + ,REP*1,FIC(6)*1
+
+ INTEGER CHOIX,I
+c + ,CONDAM
+c + ,CONDAV
+c + ,REGIME
+ LOGICAL TESTCH
+ INTEGER NBCOUP(NCOURB)
+c INTEGER OPTFPC
+c DOUBLE PRECISION ROS,DIAM,ETEN
+
+ COMMON/ETUDE/ETUDE
+ COMMON/CONST/IDPARS,IE
+ COMMON/EPSI/EPS
+ COMMON/NBVALE/NBCOUP
+
+ EPS=0.001
+ WRITE(*,*) ' '
+ WRITE(*,*) ' '
+ WRITE(*,*) ' '
+ WRITE(*,*) ' ******************************************'
+ WRITE(*,*) ' PROGRAMME DE SAISIE ET DE VERIFICATION'
+ WRITE(*,*) ' DES DONNEES ASSOCIEES AU PROGRAMME'
+ WRITE(*,*) ' DE SIMULATION DE RUPTURE DE BARRAGE'
+cliq WRITE(*,*)' RUBAR 3 '
+ WRITE(*,*)' RubarBE '
+ WRITE(*,*)' version du 9 NOVEMBRE 2016 '
+ WRITE(*,*) ' ******************************************'
+ WRITE(*,*) ' '
+ WRITE(*,*) ' '
+
+C 1)--> DESIGNATION DE L'ETUDE
+C------------------------------
+
+ WRITE(*,412)
+412 FORMAT(1H1)
+ WRITE(*,*)'DONNEZ EN MOINS DE 20 CARACTERES LE NOM DE L ETUDE'
+ READ(*,9998) ETUDE
+9998 FORMAT(A20)
+
+
+C 2)--> CHOIX DU TYPE DE SAISIE.IL EST A NOTER QUE LE PREMIER D'
+C ENTRE EUX DOIT ETRE IMPERATIVEMENT LE N° 1
+C--------------------------------------------------------------
+
+
+30 WRITE(*,*) ' '
+ WRITE(*,*) '---------------------------------------------'
+ WRITE(*,*) ' LE NOM DE L ETUDE EST : ',ETUDE
+ WRITE(*,*) 'le maillage doit exister si choix autre que 1 ou 2'
+ WRITE(*,*) '---------------------------------------------'
+ WRITE(*,*) ' MENU D APPEL DES PROGRAMMES DE SAISIE '
+ WRITE(*,*) '---------------------------------------------'
+ WRITE(*,*) 'CHOIX 1 : SAISIE DU MAILLAGE '
+ WRITE(*,*) 'CHOIX 2 : SAISIE DE LA GEOMETRIE '
+ WRITE(*,*) 'CHOIX 3 : SAISIE DU FROTTEMENT '
+ WRITE(*,*) 'CHOIX 4 : SAISIE DES DONNEES INITIALES '
+ WRITE(*,*) 'CHOIX 5 : SAISIE DES CONDITIONS AUX LIMITES '
+ WRITE(*,*) 'CHOIX 6 : SAISIE DES DEVERSEMENTS/APPORTS LATERAUX'
+ WRITE(*,*) 'CHOIX 7 : SAISIE DES OUVRAGES '
+ WRITE(*,*) 'CHOIX 8 : SAISIE DES PARAMETRES DE SAUVEGARDE'
+ WRITE(*,*) 'CHOIX 9 : SAISIE DES PARAMETRES DE SIMULATION'
+ WRITE(*,*) 'CHOIX 10 : fusion plusieurs monobiefs
+ :','en un projet multibief'
+ WRITE(*,*) 'CHOIX 11 : scission un projet monobief
+ :','en un projet multibief'
+ WRITE(*,*) 'CHOIX 12 : Transformation projet multibief
+ :','en un projet monobief pour visualisation resultats'
+ WRITE(*,*) 'CHOIX 0 : RETOUR AU SYSTEME '
+ WRITE(*,*) '---------------------------------------------'
+ WRITE(*,*) ' '
+ WRITE(*,*) 'DONNEZ VOTRE CHOIX'
+ READ (*,*) CHOIX
+
+C -- TEST DU CHOIX
+
+ TESTCH=.FALSE.
+ DO 40 I=0,12
+ TESTCH=TESTCH.OR.(CHOIX.EQ.I)
+40 CONTINUE
+ IF (.NOT.TESTCH) THEN
+ WRITE(*,*) 'VOUS AVEZ DONNE UN MAUVAIS CHOIX'
+ WRITE(*,*) 'RECOMMENCEZ'
+ GOTO 30
+ ELSEIF(CHOIX.ne.0)THEN
+
+ GOTO(1,2,3,4,5,6,7,8,9,10,11,12),CHOIX
+
+1 CALL MAILLA
+ GO TO 30
+
+2 CALL GEOMET
+ GO TO 30
+
+3 CALL FROTTE
+ GO TO 30
+
+4 CALL DONINI
+ GO TO 30
+
+5 CALL CONDLI
+ GO TO 30
+
+6 CALL SAIDEV
+ GO TO 30
+
+7 CALL SAIOUV
+ GOTO 30
+
+8 CALL SAISAV
+ GOTO 30
+
+C 3)--> REMPLISSAGE DU FICHIER DONNEE.ETUDE
+C-----------------------------------------------------------------C
+C Il contient tous les parametres de simulation C
+C-----------------------------------------------------------------C
+
+
+ 9 CALL SAICONFIG
+ GOTO 30
+
+ 10 CALL FU1BMB
+ GOTO 30
+
+ 11 CALL TR1BMB
+ GOTO 30
+
+ 12 CALL TRMB1B
+ GOTO 30
+
+
+C fin du if sur choix
+ ENDIF
+
+ END
+
+C-----------------------------------------------------------------------
+ SUBROUTINE FU1BMB
+C-----------------------------------------------------------------------
+C rassembler des projets monobiefs en un projet multibief'
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX,NBMAX,CSMAX,LNCMAX
+ PARAMETER (LMAX=3000,NBMAX=150,CSMAX=10,LNCMAX=130000)
+ CHARACTER ETUDE*20,NOMFIC*38,CHAINE*50
+ DOUBLE PRECISION X(LMAX)
+
+ COMMON/ETUDE/ETUDE
+ Write (*,*)'non operationnel'
+ return
+ end
+C-----------------------------------------------------------------------
+ SUBROUTINE TR1BMB
+C-----------------------------------------------------------------------
+C sectionner un projet monobief en un projet multibief'
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX,NBMAX,CSMAX,LNCMAX
+ PARAMETER (LMAX=3000,NBMAX=150,CSMAX=10,LNCMAX=130000)
+ CHARACTER ETUDE*20,NOMFIC*38,CHAINE*50
+ DOUBLE PRECISION X(LMAX)
+
+ COMMON/ETUDE/ETUDE
+
+ Write (*,*)'non operationnel'
+ return
+ end
+C-----------------------------------------------------------------------
+ SUBROUTINE TRMB1B
+C-----------------------------------------------------------------------
+C modifier a minima projet multibief en projet monobief
+C sert uniquement a la visualisation des resultats par pamhyr
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX,NBMAX,CSMAX,LNCMAX
+ PARAMETER (LMAX=3000,NBMAX=150,CSMAX=10,LNCMAX=130000)
+ CHARACTER ETUDE*20,NOMFIC*38,CHAINE*50
+ CHARACTER LIGNE*630,FINLIGNE*600,MORCLIGNE*58
+ DOUBLE PRECISION X(LMAX),X2(LMAX),XX
+ INTEGER NBB,M,ILM(0:NBMAX),I,NM,XNC(0:LMAX)
+ :,NSM,II,J,K,XNBCS(LNCMAX),NUMCS,NBSECT,N
+ DOUBLE PRECISION XYCOU(LNCMAX),XZCS(LNCMAX,CSMAX)
+ : ,XDCS(LNCMAX,CSmax),XSCS(LNCMAX,CSMAX),XTMCS(LNCMAX,CSMAX)
+ CHARACTER*1 XCCOU(LNCMAX)
+
+ COMMON/ETUDE/ETUDE
+
+C reecriture de geomac-i et mail
+C on lit puis on sauve avant de reecrire en ecrasant
+
+C Lecture du fichier de maillage
+C-----------------------------------------------------------------------
+ NBB=1
+ ILM(0)=0
+ OPEN(70,FILE = 'confl.'//ETUDE,STATUS='OLD',ERR=111)
+ READ(70,*,ERR=111)NBB
+ CLOSE(70)
+ NOMFIC='mail.'//ETUDE
+ OPEN(65,FILE=NOMFIC,STATUS='UNKNOWN')
+ DO 1000 M=1,NBB
+ READ(65,*) ILM(M)
+ ILM(M)=ILM(M-1)+ILM(M)
+ READ(65,*) (X(I),I=ILM(M-1)+1,ILM(M))
+ READ(65,*) (X2(I),I=ILM(M-1)+1,ILM(M)-1)
+ X2(ILM(m))=X2(ILM(m)-1)
+ 1000 CONTINUE
+ CLOSE(65)
+ NM=ILM(NBB)
+ GO TO 112
+ 111 NBB=1
+ Write(*,*)'un seul bief initialement'
+ write(*,*)'erreur'
+ return
+C Ouverture du fichier de géométrie abscisse-cote aux intermailles
+C geomac-i
+ 112 NOMFIC='geomac-i.'//ETUDE
+ OPEN(24,FILE=NOMFIC,STATUS='OLD',FORM='FORMATTED',ERR=98)
+ XNC(0)=0
+ GO TO 980
+ 98 WRITE(*,*)'FICHIER geomac-i INEXISTANT'
+ Stop
+ 980 DO M=1,NBB
+ read(24,*)nsm
+ IF(NSM.NE.ILM(M)-1-ILM(M-1))THEN
+ write(*,*)'erreur en tete bief geomac-i ou fichier mail'
+ write(*,*)'nombre sections bief = ',m
+ Stop
+ ENDIF
+
+ DO I=ILM(M-1)+1,ILM(M)-1
+C Contrôle des en-têtes de sections supprime
+ READ(24,*,END=240) II,XX,XNC(I)
+ XNC(I)=XNC(I)+XNC(I-1)
+ DO J=XNC(I-1)+1,XNC(I)
+C Lecture des couples abscisse-cote
+ READ(24,'(A630)') LIGNE
+ READ(LIGNE,9994) XCCOU(J),XYCOU(J),FINLIGNE
+C Décodage sédimentaire (description point par point)
+ DO NUMCS=1,CSMAX
+ MORCLIGNE=FINLIGNE(1+58*(NUMCS-1):58*NUMCS)
+ READ(MORCLIGNE,9996,END=200) XZCS(J,NUMCS),XDCS(J,NUMCS)
+ & ,XSCS(J,NUMCS),XTMCS(J,NUMCS)
+ ENDDO
+ 200 CONTINUE
+ XNBCS(J)=NUMCS-1
+C fin boucle sur J
+ ENDDO
+C fin boucle sur I
+ ENDDO
+ XNC(ILM(M))=XNC(ILM(M)-1)
+C fin boucle sur IB
+
+ ENDDO
+ 240 CLOSE(24)
+
+ 9994 FORMAT(A1,1X,F11.5,A600)
+ 9996 FORMAT(F13.5,2F15.10,F15.5)
+
+
+ NOMFIC='mail.'//ETUDE//'ancien'
+ OPEN(65,FILE=NOMFIC,STATUS='UNKNOWN')
+ DO 1015 M=1,NBB
+ write(65,*) ILM(M)-ILM(m-1)
+ Write(65,*) (X(I),I=ILM(M-1)+1,ILM(M))
+ WRITE(65,*) (X2(I),I=ILM(M-1)+1,ILM(M)-1)
+ 1015 CONTINUE
+ CLOSE(65)
+
+ NOMFIC='mail.'//ETUDE
+ OPEN(65,FILE=NOMFIC,STATUS='UNKNOWN')
+ write(65,*) ILM(M)
+ DO 1001 M=1,NBB
+ Write(65,*) (X(I),I=ILM(M-1)+1,ILM(M))
+ 1001 CONTINUE
+ DO 1002 M=1,NBB-1
+ WRITE(65,*) (X2(I),I=ILM(M-1)+1,ILM(M)-1)
+ WRITE(65,*) X2(ILM(M)-1)
+ 1002 CONTINUE
+ M=NBB
+ WRITE(65,*) (X2(I),I=ILM(M-1)+1,ILM(M)-1)
+ CLOSE(65)
+
+C RéÉcriture de 'geomac-i.etude' initial
+C-----------------------------------------------------------------------
+ NOMFIC='geomac-i.'//ETUDE//'ancien'
+ OPEN(24,FILE=NOMFIC,STATUS='UNKNOWN')
+ DO M=1,NBB
+ NBSECT=ILM(M)-ILM(m-1)-1
+ WRITE(24,110)NBSECT,0.
+ DO I=ilm(m-1)+1,ilm(m-1)+NBSECT
+ WRITE(24,2111)I,X2(I),xnc(i)-xnc(i-1)
+ DO J=XNC(I-1)+1,XNC(I)
+ DO K=1,XNBCS(J)
+ N=1+58*(K-1)
+ WRITE(FINLIGNE(N:),113)xZCS(J,K),xDCS(J,K),xSCS(J,K)
+ & ,xTMCS(J,K)
+ ENDDO
+ WRITE(24,2112)xcCOU(J),xYCOU(J),FINLIGNE
+ ENDDO
+ ENDDO
+C fin boucle sur bief
+ ENDDO
+ CLOSE(24)
+ WRITE(*,*)'Ecriture de ',NOMFIC
+
+C Écriture de 'geomac-i.etude'
+C-----------------------------------------------------------------------
+ NOMFIC='geomac-i.'//ETUDE
+ OPEN(24,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(24,110)ilm(nbb)-1,0.
+ DO M=1,NBB
+ NBSECT=ILM(M)-ILM(m-1)-1
+ DO I=ilm(m-1)+1,ilm(m-1)+NBSECT
+ WRITE(24,2111)I,X2(I),xnc(i)-xnc(i-1)
+ DO J=XNC(I-1)+1,XNC(I)
+ DO K=1,XNBCS(J)
+ N=1+58*(K-1)
+ WRITE(FINLIGNE(N:),113)xZCS(J,K),xDCS(J,K),xSCS(J,K)
+ & ,xTMCS(J,K)
+ ENDDO
+ WRITE(24,2112)xcCOU(J),xYCOU(J),FINLIGNE
+ ENDDO
+ ENDDO
+ IF(M.NE.NBB)THEN
+C on double la dernieres section de chaque bief sauf le dernier
+ I=ilm(m-1)+NBSECT
+ WRITE(24,2111)I,X2(I),xnc(i)-xnc(i-1)
+ DO J=XNC(I-1)+1,XNC(I)
+ DO K=1,XNBCS(J)
+ N=1+58*(K-1)
+ WRITE(FINLIGNE(N:),113)xZCS(J,K),xDCS(J,K),xSCS(J,K)
+ & ,xTMCS(J,K)
+ ENDDO
+ WRITE(24,2112)xcCOU(J),xYCOU(J),FINLIGNE
+ ENDDO
+C fin du endif sur m=NBB
+ ENDIF
+C fin boucle sur bief
+ ENDDO
+ CLOSE(24)
+ WRITE(*,*)'Ecriture de ',NOMFIC
+ 110 FORMAT(I4,1X,F12.3)
+ 2111 FORMAT(I4,1X,F11.3,1X,I4)
+ 2112 FORMAT(A1,1X,F11.5,A600)
+ 113 FORMAT(F13.5,F15.10,F15.10,F15.5)
+C 2111 FORMAT(3(A7,F16.3))
+C 2112 FORMAT(1X,I4,1X,3(F16.3))
+
+ RETURN
+ END
+
+C-----------------------------------------------------------------------
+ SUBROUTINE SAIDONNEE
+C-----------------------------------------------------------------------
+C Remplissage du fichier des parametres de simulation 'donnee.etude'
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER CONDAM,CONDAV,REGIME,OPTFPC
+ CHARACTER ETUDE*20,NOMFIC*38,CHAINE*50
+ & ,IODEV*1,IOPDT*1,REP*1,IOVIS*1,FIC(7)*1
+ DOUBLE PRECISION CFL,TINIT,TMAX
+ & ,TIOPDT,DT
+ & ,TS,PSAVE,DTSAUV
+ & ,FDEB1,FDEB2,FDEB3
+ :,ROS,DIAM,ETEN
+
+ COMMON/ETUDE/ETUDE
+
+C Récupération des informations
+C-----------------------------------------------------------------------
+ NOMFIC='donnee.'//ETUDE
+ OPEN(10,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(*,*) 'REMPLISSAGE DU FICHIER DONNEE.ETUDE CONTENANT'
+ WRITE(*,*) ' LES PARAMETRES PERMETTANT '
+ WRITE(*,*) ' D ''INITIALISER LES CALCULS '
+ WRITE(*,*)
+ WRITE(*,*)
+ WRITE(*,*) 'RENTREZ LES PARAMETRES DE CALCUL :'
+ WRITE(*,*)
+ WRITE(*,*) ' CONDITION DE FRIEDRICH LEVY (CFL)'
+ READ(*,*) CFL
+ CHAINE='CONDITION DE FRIEDRICH LEVY='
+ WRITE(10,101)CHAINE,CFL
+ WRITE(*,*)'- Condition limite amont'
+ WRITE(*,*)' - Rien ---------------> 0'
+ WRITE(*,*)' - Q(t)[hydro] --------> 1'
+ WRITE(*,*)' - Z(t)[hamont] -------> 2'
+ READ(*,*) CONDAM
+ CHAINE='CONDITION LIMITE AMONT='
+ WRITE(10,102)CHAINE,CONDAM
+ WRITE(*,*) ' CONDITION LIMITE AVAL '
+ WRITE(*,*) ' REGIME CRITIQUE-----> 7 '
+ WRITE(*,*) ' REGIME UNIFORME-----> 5 '
+ WRITE(*,*) ' Z(T) ---- ----------> 3 '
+ WRITE(*,*) ' REFLEXION ----------> 2 '
+ WRITE(*,*) ' Q(Y) ---------------> 1 '
+ WRITE(*,*) ' SORTIE LIBRE -------> 0 '
+ READ(*,*) CONDAV
+ CHAINE='CONDITION LIMITE AVAL='
+ WRITE(10,103)CHAINE,CONDAV
+103 FORMAT(A50,I1)
+ WRITE(*,*) ' COTE IMPOSEE A L''AMONT '
+ WRITE(*,*) ' NON ---->0 '
+ WRITE(*,*) ' OUI ---->1 '
+C WRITE(*,*) ' OUI SEULE ---->2 '
+ READ(*,*)REGIME
+ CHAINE='COTE IMPOSEE EN AMONT='
+ WRITE(10,119)CHAINE,REGIME
+119 FORMAT(A50,I1)
+ WRITE(*,414)
+414 FORMAT(1H1)
+ WRITE(*,*) 'RENTREZ LES OPTIONS DE CALCUL : '
+ WRITE(*,*)
+ WRITE(*,*) ' DEVERSEMENT LATERAL '
+ WRITE(*,*) ' NON ----> N OUI ----> O '
+ READ(*,'(A1)') IODEV
+ CHAINE='DEVERSEMENT LATERAL O/N='
+ WRITE(10,104)CHAINE,IODEV
+ WRITE(*,*) ' PAS DE TEMPS AUTOMATIQUE '
+ WRITE(*,*) ' NON -->N OUI --> O MIXTE --> M '
+ READ(*,'(A1)') IOPDT
+ CHAINE='PAS DE TEMPS AUTOMATIQUE O/N/M'
+ WRITE(10,105)CHAINE,IOPDT
+105 FORMAT(A50,A1)
+ WRITE(*,*) ' CREATION DU FICHIER largeur et autres '
+ WRITE(*,*) ' NON -->N OUI --> O '
+ READ(*,'(A1)') IOVIS
+ CHAINE='CREATION FICHIER largeur O/N='
+ WRITE(10,106)CHAINE,IOVIS
+106 FORMAT(A50,A1)
+ WRITE(*,415)
+415 FORMAT(1H1)
+ WRITE(*,*) 'RENTREZ LES PARAMETRES TEMPORELS '
+ WRITE(*,*)
+ WRITE(*,*) ' REPRISE DU CALCUL '
+ WRITE(*,*) ' NON ----> N OUI ----> O '
+ READ(*,'(A1)') REP
+ CHAINE='REPRISE DU CALCUL A T DIFFERENT DE 0'
+ WRITE(10,107)CHAINE,REP
+107 FORMAT(A50,A1)
+ WRITE(*,*) 'TEMPS DEBUT DE CALCUL '
+ READ(*,*) TINIT
+ CHAINE='TEMPS DEBUT DE CALCUL='
+ WRITE(10,108)CHAINE,TINIT
+108 FORMAT(A50,E12.5)
+ WRITE(*,*) ' TEMPS DE FIN CALCUL '
+ READ(*,*) TMAX
+ CHAINE='TEMPS FIN DE CALCUL='
+ WRITE(10,109)CHAINE,TMAX
+109 FORMAT(A50,E12.5)
+ IF(IOPDT.NE.'M')THEN
+ WRITE(*,*) ' TEMPS JUSQU''AUQUEL DT=CONSTANT '
+C WRITE(*,*) ' METTEZ 0 SI PAS DE TEMPS MIXTE PAS CHOISI'
+ READ(*,*) TIOPDT
+ ELSE
+ TIOPDT=TINIT
+ ENDIF
+ CHAINE='PAS DE TEMPS CONSTANT JUSQU''A T='
+ WRITE(10,110)CHAINE,TIOPDT
+110 FORMAT(A50,E12.5)
+ WRITE(*,*) ' PAS DE TEMPS INITIAL DT '
+ READ(*,*) DT
+ CHAINE='VALEUR DU PAS DE TEMPS INITIAL ='
+ WRITE(10,111)CHAINE,DT
+111 FORMAT(A50,E12.5)
+ WRITE(*,416)
+416 FORMAT(1H1)
+ WRITE(*,417)
+417 FORMAT(1H1)
+ WRITE(*,*) 'FICHIERS RESULTATS DESIRES '
+ WRITE(*,*) ' NON ----> N OUI ----> O '
+ WRITE(*,*)
+ WRITE(*,*) ' tps (REPRISE ULTERIEURE) '
+ READ(*,'(A1)') FIC(1)
+ WRITE(*,*) ' hydlim (SAUVEGARDE Q,V,H à x fixé )'
+ READ(*,'(A1)') FIC(2)
+ WRITE(*,*) ' envlop (ENVELOPPE DES Q,Y,V) '
+ READ(*,'(A1)') FIC(3)
+ WRITE(*,*) ' profil (SAUVEGARDE A T DONNE) '
+ READ(*,'(A1)') FIC(4)
+ WRITE(*,*) ' profil (SAUVEGARDE dt donne) '
+ READ(*,'(A1)') FIC(5)
+ WRITE(*,*) ' trajec (TEMPS ARRIVEE DES FRONTS) '
+ READ(*,'(A1)') FIC(6)
+ WRITE(*,*) 'RENTREZ LES PARAMETRES DE SORTIE '
+ WRITE(*,*)
+ If(fic(1).eq.'O'.or.fic(1).eq.'o')then
+ WRITE(*,*) ' TEMPS DE SAUVEGARDE POUR REPRISE EVENTUELLE '
+ READ(*,*) TS
+ else
+ ts=tmax
+ endif
+ CHAINE='TEMPS DE SAUVEGARDE POUR REPRISE='
+ WRITE(10,112)CHAINE,TS
+112 FORMAT(A50,E12.5)
+ WRITE(*,*)
+ +' PAS DE SAUVEGARDE POUR LES SORTIES DE LIGNES D''EAU'
+ READ(*,*) DTSAUV
+ CHAINE='PAS DE SAUVEGARDE POUR LIGNE D''EAU='
+ WRITE(10,113)CHAINE,DTSAUV
+ if(fic(2).eq.'O'.or.fic(2).eq.'o')then
+ WRITE(*,*)
+ +' PAS DE SAUVEGARDE POUR LES SORTIES D''HYDROGRAMMES'
+ READ(*,*) PSAVE
+ else
+ psave=tmax
+ endif
+ CHAINE='PAS DE SAUVEGARDE POUR HYDROGRAMMES'
+ WRITE(10,113)CHAINE,PSAVE
+113 FORMAT(A50,E12.5)
+ if(fic(6).eq.'O'.or.fic(6).eq.'o')then
+ WRITE(*,*) ' '
+ WRITE(*,*) ' ENTREZ LES TROIS FRONTS DE DEBIT '
+ WRITE(*,*) ' FRONT DEBIT 1 '
+ WRITE(*,*) ' FRONT DEBIT 2 '
+ WRITE(*,*) ' FRONT DEBIT 3 '
+ READ(*,*) FDEB1,FDEB2,FDEB3
+ else
+ fdeb1=1.
+ fdeb2=10.
+ fdeb3=100.
+ endif
+ CHAINE='FRONTS DE DEBIT='
+ WRITE(10,114)CHAINE,FDEB1
+ WRITE(10,114)CHAINE,FDEB2
+ WRITE(10,114)CHAINE,FDEB3
+114 FORMAT(A50,E12.5)
+ CHAINE='tps (REPRISE)-------------> O/N'
+ WRITE(10,115)CHAINE,FIC(1)
+ CHAINE='hydlim (Q(T),V(T),H(T))---> O/N'
+ WRITE(10,115)CHAINE,FIC(2)
+ CHAINE='envlop (QMAX,YMAX,VMAX)---> O/N'
+ WRITE(10,115)CHAINE,FIC(3)
+ CHAINE='profil(Y(X),V(X),Q(X) A T)--> O/N'
+ WRITE(10,115)CHAINE,FIC(4)
+ CHAINE='profil(Y(X),V(X),Q(X) A DT)--> O/N'
+ WRITE(10,115)CHAINE,FIC(5)
+ CHAINE='trajec (FRONTS)-----------> O/N'
+ WRITE(10,115)CHAINE,FIC(6)
+ WRITE(*,*) ' transport de sediments O/N'
+ READ(*,'(A1)') FIC(6)
+ CHAINE='transport de sediment (o,n)'
+cliq FIC(6)='N'
+ WRITE(10,115)CHAINE,FIC(6)
+
+115 FORMAT(A50,A1)
+ WRITE(*,*) ' OPTION FROTTEMENTS '
+ WRITE(*,*) ' pas de changement --> 0 '
+ WRITE(*,*) ' BROWNLIE -----> 1 '
+ WRITE(*,*) ' GRIFFITH ----------> 2 '
+ WRITE(*,*) ' KARIM ----------> 3 '
+ WRITE(*,*) ' WU et WANG ---------> 4 '
+ WRITE(*,*) ' YU et LIM -------> 5 '
+ WRITE(*,*) ' RECKING -------> 6 '
+ WRITE(*,*) ' GARDE et RAJU ------> 7 '
+ READ(*,*) OPTFPC
+ CHAINE='OPTION FROTTEMENTS='
+ WRITE(10,103)CHAINE,OPTFPC
+ IF(OPTFPC.NE.0)THEN
+ WRITE(*,*) ' masse volumique des sediments '
+ READ(*,*) ROS
+ CHAINE='MASSE VOLUMIQUE DES SEDIMENTS='
+ WRITE(10,109)CHAINE,ROS
+ WRITE(*,*) ' diametre des sediments '
+ READ(*,*) DIAM
+ CHAINE='DIAMETRE DES SEDIMENTS='
+ WRITE(10,109)CHAINE,DIAM
+ WRITE(*,*) ' etendue des sediments '
+ READ(*,*) ETEN
+ CHAINE='ETENDUE DES SEDIMENTS='
+ WRITE(10,109)CHAINE,ETEN
+C si optfpc =0
+ ELSE
+ CHAINE='MASSE VOLUMIQUE DES SEDIMENTS='
+ WRITE(10,109)CHAINE,2600.
+ CHAINE='DIAMETRE DES SEDIMENTS='
+ WRITE(10,109)CHAINE,0.1
+ CHAINE='ETENDUE DES SEDIMENTS='
+ WRITE(10,109)CHAINE,1.
+ ENDIF
+ CLOSE(10)
+
+ RETURN
+101 FORMAT(A50,E12.5)
+102 FORMAT(A50,I1)
+104 FORMAT(A50,A1)
+
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE SAIVISU
+C-----------------------------------------------------------------------
+C Remplissage du fichier des parametres de simulation 'config-visu'
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ LOGICAL TRGEOMD,TRGEOMI,TRGEOML,TRSURFL,TRSL,ASPSL
+ CHARACTER ETUDE*20,NOMFIC*38,REP*1,CHAINE*50,COMM*60
+ DOUBLE PRECISION XMIN,XMAX,DIRVUE,XVUE,YVUE,ZVUE,LARVUE
+
+ COMMON/ETUDE/ETUDE
+
+C Récupération des informations
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' REMPLISSAGE DU FICHIER config-visu '
+ WRITE(*,*)' DES PARAMETRES DE VISUALISATION 3D EN VRML '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)'Le fichier de sortie ''visu.wrl'' represente une vue '
+ & ,'3D visitable du cas modelise. Il est lisible par tout '
+ & ,'navigateur muni d''un plug-in gratuit de gestion du VRML.'
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ 100 CONTINUE
+ WRITE(*,*)'Quel est le minimum des abscisses en long (m) ?'
+ READ(*,*)XMIN
+ WRITE(*,*)'Quel est le maximum des abscisses en long (m) ?'
+ READ(*,*)XMAX
+ IF(XMAX.LE.XMIN)THEN
+ WRITE(*,*)'PROBLEME d''ordre des abscisses. Recommencez.'
+ GOTO 100
+ ENDIF
+ WRITE(*,*)'Quelles sont les coordonnées X,Y,Z du point '
+ & ,'d''observation dans le repere des points du profil (m m m) ?'
+ READ(*,*)XVUE,YVUE,ZVUE
+ WRITE(*,*)'Quelle est la direction d''observation (radians) ?'
+ WRITE(*,*)'[axeX:-1.57 axe-Y:0.00 axe-X:1.57 axeY:3.14]'
+ READ(*,*)DIRVUE
+ WRITE(*,*)'Quelle est la largeur du champ de vision (radians) ? '
+ & ,'[0.35]'
+ READ(*,*)LARVUE
+ WRITE(*,*)'Voulez-vous tracer la topographie du lit '
+ & ,'(strate superieure) ? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ TRGEOML=.TRUE.
+ WRITE(*,*)'Voulez-vous materialiser sa surface ? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ TRSURFL=.TRUE.
+ ELSE
+ TRSURFL=.FALSE.
+ ENDIF
+ ELSE
+ TRGEOML=.FALSE.
+ TRSURFL=.FALSE.
+ ENDIF
+ WRITE(*,*)'Voulez-vous tracer la topographie de la strate '
+ & ,'inferieure du lit (eventuellement inerodable) ? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ TRGEOMD=.TRUE.
+ ELSE
+ TRGEOMD=.FALSE.
+ ENDIF
+ WRITE(*,*)'Voulez-vous tracer la topographie du lit initial '
+ & ,'? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ TRGEOMI=.TRUE.
+ ELSE
+ TRGEOMI=.FALSE.
+ ENDIF
+ WRITE(*,*)'Voulez-vous tracer la surface libre ? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ TRSL=.TRUE.
+ WRITE(*,*)'Voulez-vous moirer sa surface ? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ ASPSL=.TRUE.
+ ELSE
+ ASPSL=.FALSE.
+ ENDIF
+ ELSE
+ TRSL=.FALSE.
+ ASPSL=.FALSE.
+ ENDIF
+ WRITE(*,*)'Quel est le titre ? (60 caracteres)'
+ READ(*,*)COMM
+
+C Écriture dans 'config-visu'
+C-----------------------------------------------------------------------
+ NOMFIC='config-visu'
+ OPEN(10,FILE=NOMFIC,STATUS='UNKNOWN')
+ CHAINE='Minimum des abscisses en long (m)'
+ WRITE(10,101)XMIN,CHAINE
+ CHAINE='Maximum des abscisses en long (m)'
+ WRITE(10,101)XMAX,CHAINE
+ CHAINE='Compartiment inférieur courant'
+ IF(TRGEOMD)THEN
+ WRITE(10,102)1,CHAINE
+ ELSE
+ WRITE(10,102)0,CHAINE
+ ENDIF
+ CHAINE='Compartiment supérieur initial'
+ IF(TRGEOMI)THEN
+ WRITE(10,102)1,CHAINE
+ ELSE
+ WRITE(10,102)0,CHAINE
+ ENDIF
+ CHAINE='Compartiment supérieur courant'
+ IF(TRGEOML)THEN
+ WRITE(10,102)1,CHAINE
+ ELSE
+ WRITE(10,102)0,CHAINE
+ ENDIF
+ CHAINE='Surface du compartiment supérieur courant'
+ IF(TRSURFL)THEN
+ WRITE(10,102)1,CHAINE
+ ELSE
+ WRITE(10,102)0,CHAINE
+ ENDIF
+ CHAINE='Surface libre / aspect moiré'
+ IF(TRSL.AND.ASPSL)THEN
+ WRITE(10,103)1,1,CHAINE
+ ELSEIF(TRSL.AND..NOT.ASPSL)THEN
+ WRITE(10,103)1,0,CHAINE
+ ELSE
+ WRITE(10,103)0,0,CHAINE
+ ENDIF
+ CHAINE='Dir. obs. (rad) [X:-1.57 -Y:0.00 -X:1.57 Y:3.14]'
+ WRITE(10,101)DIRVUE,CHAINE
+ CHAINE='X,Y,Z observation (m m m)'
+ WRITE(10,104)XVUE,YVUE,ZVUE,CHAINE
+ CHAINE='Largeur du champ de vision (rad)'
+ WRITE(10,101)LARVUE,CHAINE
+ WRITE(10,105)COMM
+ CLOSE(10)
+
+ RETURN
+101 FORMAT(F11.3,1X,A50)
+102 FORMAT(3X,I1,8X,A50)
+103 FORMAT(3X,I1,3X,I1,4X,A50)
+104 FORMAT(3(F11.3,1X),A50)
+105 FORMAT(A60)
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE SAITS
+C-----------------------------------------------------------------------
+C Remplissage du fichier des parametres du transport solide 'ts.etude'
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER OPTS,ODCHAR,UNISOL,TYPDEF,CHOIXC,OPTION
+ & ,DEPOT,CAPASOL,DEFOND,VARCONS
+ INTEGER DEMIX
+ CHARACTER ETUDE*20,NOMFIC*38,CHAINE*50
+ & ,VERSION*6
+ DOUBLE PRECISION ROS,POR,VISC
+ & ,DCHARG,HALFA,BMIU,DCHARD,DCHARS
+ & ,MUCASO,TCADIM
+
+ COMMON/ETUDE/ETUDE
+
+C Récupération des informations
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' REMPLISSAGE DU FICHIER ts.etude CONTENANT '
+ WRITE(*,*)' LES PARAMETRES DU TRANSPORT SOLIDE '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'Rentrez les options de calcul:'
+C Le menu de saisie des conditions initiales solides ne permet pas
+C d'entrer des concentrations. Donc on reste homogène dans cette version
+ WRITE(*,*)'- Unites pour les solides'
+ WRITE(*,*)' - Qs en kg/s ---------> 1'
+ WRITE(*,*)' - Cs en kg/m3 --------> 2'
+ READ(*,*)UNISOL
+C UNISOL=1
+ WRITE(*,*)'- Type de loi de capacite solide'
+ WRITE(*,*)' - Meyer-Peter&Muller -> 1'
+ WRITE(*,*)' - Engelund&Hansen ----> 2'
+ WRITE(*,*)' - Bagnold ------------> 3'
+ WRITE(*,*)' - Meyer-Peter&Muller+correction contrainte efficace '
+ & ,'----------> 4'
+ WRITE(*,*)' - Meyer-Peter&Muller+contrainte critique '
+ & ,'adimensionnelle fixe ----> 5'
+ WRITE(*,*)' - Meyer-Peter&Muller+correction contrainte efficace+'
+ WRITE(*,*)' contrainte critique adimensionnelle '
+ & ,'fixe -----------> 6'
+ WRITE(*,*)' - Meyer-Peter&Muller modifie pour cours d eau '
+ & ,'polonais ------------> 7'
+ WRITE(*,*)' - Ackers White -> 8'
+ WRITE(*,*)' - Smart Jaeggi -> 9'
+ WRITE(*,*)' - van Rijn -> 10'
+ WRITE(*,*)' - Rickenmann -> 11'
+ WRITE(*,*)' - Camenen et Larson -> 12'
+ WRITE(*,*)' - Schoklitsch -> 13'
+ WRITE(*,*)' - Sato et al. -> 14'
+ WRITE(*,*)' - Recking et al. -> 15'
+ READ(*,*)OPTS
+ WRITE(*,*)'- Option pour la distance de chargement'
+ WRITE(*,*)' - Fixe ---------------> 0'
+ WRITE(*,*)' - Calcul loi de Han --> 1'
+ write(*,*)' - Wu et Wang ---> 2'
+ READ(*,*)ODCHAR
+C Pour l'instant, la méthode de déformation par les contraintes locales
+C (Méthode des Perpendiculaires Confondues) n'est pas disponible
+C WRITE(*,*)'- Methode de deformation geometrique de la section'
+C WRITE(*,*)' - Synthetique --------> 1'
+C WRITE(*,*)' - Contrainte locale --> 2'
+C READ(*,*) TYPDEF
+ TYPDEF=1
+ WRITE(*,*)''
+ WRITE(*,*)' en cas de depot,vous pouvez choisir'
+ &,'trois facons de l appliquer'
+ WRITE(*,*)''
+ WRITE(*,*)'- DEPOT 1: dans le lit mineur par couches '
+ & ,'horizontales------> 1'
+ WRITE(*,*)''
+ WRITE(*,*)'- DEPOT 2 :dans le lit mineur, depot '
+ & ,'uniforme ---------> 2'
+ WRITE(*,*)''
+ WRITE(*,*)'- DEPOT 3: dans le lit mineur, proportionnel'
+ & ,' a 1/contrainte ----------->3'
+ WRITE(*,*)''
+ WRITE(*,*)'- DEPOT 4: dans le lit mineur,proportionnel a '
+ & ,'(contrainte critique-K*contrainte)^(m)'
+ WRITE(*,*)' K est fonction de min(contrainte critique'
+ & ,'/contrainte) -------------->4'
+ WRITE(*,*)''
+ WRITE(*,*)'- DEPOT 5: dans le lit mineur,proportionnel a '
+ & ,'(contrainte critique-K*contrainte)^(m)'
+ WRITE(*,*)' K est fonction de moyenne(contrainte critique'
+ & ,'/contrainte) -------------->5'
+ WRITE(*,*)''
+ WRITE(*,*)'- DEPOT 6: dans le lit mineur,proportionnel a '
+ & ,'(contrainte)^(m)'
+ WRITE(*,*)' m est fonction de la formule de capacité solide '
+ & ,'utilisee---------------->6'
+ READ(*,*) DEPOT
+ WRITE(*,*)''
+ WRITE(*,*)' Avec quelle methode vous voulez calculer '
+ & ,'la contrainte hydrodynamique:'
+ WRITE(*,*)' - Choix 1:MPC ----------------------> 1'
+ WRITE(*,*)' - Choix 2: formule regime uniforme---------> 2'
+ WRITE(*,*)' - Choix 3:MPC pente moyenne sur 3 points --> 3'
+ WRITE(*,*)' - Choix 4:MPC moyenne sur 3 segments ------> 4'
+ WRITE(*,*)' - Choix 5:MPC moyenne sur lp=h ------------> 5'
+ WRITE(*,*)' - Choix 6:MPC selon diametre coucheactive -> 6'
+ READ(*,*) CHOIXC
+ WRITE(*,*)' Avec quelle methode vous voulez calculer la '
+ & ,'contrainte critique:'
+ WRITE(*,*)' - Option 1:Contrainte critique avec facteur IKDEA '
+ & ,'-------> 1'
+ WRITE(*,*)' - Option 2: Contrainte critique sans facteur IKEDA '
+ & ,'------------------> 2'
+ READ(*,*) OPTION
+ WRITE(*,*)''
+ WRITE(*,*)' Avec quelle methode vous voulez calculer '
+ & ,'la capacite solide:'
+ WRITE(*,*)' - Choix 1:largeur active et debit solide UNITAIRE '
+ & ,'sont supposes independants de la repartition des contraintes'
+ WRITE(*,*)''
+ WRITE(*,*)'Attention:cette methode est valable uniquement si vous'
+ & ,'avez choisi de calculer la contrainte par la formule "regime'
+ & ,'uniforme" et si la contrainte critique ne depend pas'
+ & ,'de la pente ----------------------> 1'
+ WRITE(*,*)''
+ WRITE(*,*)' - Choix2: largeur active calculee en fonction de la'
+ & ,'repartition de la contrainte mais la capacite solide UNITAIRE'
+ & ,'est fonction de la contrainte moyenne---------> 2'
+ WRITE(*,*)''
+ WRITE(*,*)' - Choix3: largeur active et capacite solide UNITAIRE'
+ & ,'sont fonction de la repartition des contraintes---------> 3'
+ READ(*,*) CAPASOL
+ WRITE(*,*)''
+ WRITE(*,*)'Rentrez les parametres:'
+ WRITE(*,*)'- Viscosite cinematique de l eau (souvent 1.3E-6 m2/s)'
+ READ(*,*)VISC
+ WRITE(*,*)'CONTRAINTE CRITIQUE ADIMENSIONNELLE (souvent 0.047)'
+ read(*,*)TCADIM
+ WRITE(*,*)'- Masse volumique des sediments (kg/m3)'
+ READ(*,*)ROS
+ WRITE(*,*)'- Porosite des sediments (m3/m3)'
+ READ(*,*)POR
+ WRITE(*,*)'- tangente de l angle frottement interne des sediments'
+ READ(*,*)BMIU
+ IF(ODCHAR.EQ.0)THEN
+ WRITE(*,*)'- Distance de chargement (m)'
+ READ(*,*)DCHARG
+ HALFA=1.
+ ELSE
+ WRITE(*,*)'- Parametre de chargement de la loi de Han (m)'
+ READ(*,*)HALFA
+ DCHARG=1.
+ ENDIF
+ WRITE(*,*)' Effectuer le demixage'
+ WRITE(*,*)' - Non ----------------> 0'
+ WRITE(*,*)' - Oui ----------------> 1'
+ READ(*,*)DEMIX
+ IF(DEMIX.EQ.1)THEN
+ WRITE(*,*)'- Distance de chargement pour le demixage des '
+ :,'diametres (m)'
+ READ(*,*)DCHARD
+ WRITE(*,*)'- Distance de chargement pour le demixage des '
+ : ,'étendues granulométriques(m)'
+ READ(*,*)DCHARS
+ ENDIF
+ WRITE(*,*)'- Multiplicateur de la capacite solide (adim)'
+ READ(*,*)MUCASO
+ WRITE(*,*)' Deformation du lit'
+ WRITE(*,*)' - Non ----------------> 0'
+ WRITE(*,*)' - Oui ----------------> 1'
+ READ(*,*)DEFOND
+ WRITE(*,*)' variable hydraulique a conserver'
+ WRITE(*,*)' - section mouillee et debit------------> 1'
+ WRITE(*,*)' - niveau d eau et debit----------------> 2'
+ WRITE(*,*)' - niveau deau et vitesse--------------> 3'
+ READ(*,*)VARCONS
+C Écriture dans 'ts.etude' (au format donné par VERSION)
+C-----------------------------------------------------------------------
+c VERSION='rube1b'
+ VERSION=' '
+ NOMFIC='ts.'//ETUDE
+ OPEN(32,FILE=NOMFIC,STATUS='UNKNOWN')
+ CHAINE='Masse volumique des sediments (kg/m3)'
+ WRITE(32,101)CHAINE,ROS,VERSION
+C WRITE(32,102)CHAINE,ROS
+ CHAINE='Porosite des sediments (m3/m3)'
+ WRITE(32,102)CHAINE,POR
+ CHAINE='Distance de chargement (si fixe) (m)'
+ WRITE(32,102)CHAINE,DCHARG
+ CHAINE='Parametre de chargement de la loi de Han (m)'
+ WRITE(32,102)CHAINE,HALFA
+ CHAINE='Multiplicateur de la capacite solide (adim)'
+ WRITE(32,102)CHAINE,MUCASO
+ CHAINE='Contr.crit.adim. /viscosite '
+ WRITE(32,105)CHAINE,TCADIM,VISC
+ CHAINE='Type de loi de capacite solide (1,2,3,4,5,6,7)'
+ WRITE(32,103)CHAINE,OPTS
+ CHAINE='Option pour la distance de chargement (0,1)'
+ WRITE(32,103)CHAINE,ODCHAR
+ CHAINE='Choix debits solides / concentrations (1,2)'
+ WRITE(32,103)CHAINE,UNISOL
+ CHAINE='Precision reactualisation geometrie (0-4)'
+ WRITE(32,103)CHAINE,TYPDEF
+ CHAINE='methode dapplication du depot (1,2,3,4,5,6)'
+ WRITE(32,103)CHAINE,DEPOT
+ CHAINE='Choix de calcul de la contrainte (1,2)'
+ WRITE(32,103)CHAINE,CHOIXC
+ CHAINE='Option de calcul de la contrainte critique (1,2)'
+ WRITE(32,103)CHAINE,OPTION
+ CHAINE='methode de calcul de la capacité solide (1,2,3)'
+ WRITE(32,103)CHAINE,CAPASOL
+ CHAINE='tangente angle frottement moyen des sediments'
+ WRITE(32,102)CHAINE,BMIU
+ CHAINE='Demixage (0,1)'
+ WRITE(32,103)CHAINE,DEMIX
+ CHAINE='deformation du lit (0,1)'
+ WRITE(32,103)CHAINE,DEFOND
+ CHAINE='variables a conserver (1,2,3)'
+ WRITE(32,103)CHAINE,VARCONS
+ IF(DEMIX.NE.1)THEN
+ dchard=dcharg
+ dchars=dcharg
+ ENDIF
+ CHAINE='Distance de chargement pour les diametres de sediment (m)'
+ WRITE(32,102)CHAINE,DCHARD
+ CHAINE='Distance de chargement pour l etendue granulometrique (m)'
+ WRITE(32,102)CHAINE,DCHARS
+ CLOSE(32)
+
+
+ RETURN
+101 FORMAT(A50,F10.5,A6)
+102 FORMAT(A50,F10.5)
+103 FORMAT(A50,I1)
+104 FORMAT(A50,F10.8)
+105 format(a50,f10.5,f10.8)
+ END
+
+C-----------------------------------------------------------------------
+ SUBROUTINE SAICONFIG
+C-----------------------------------------------------------------------
+C Menu relatif à la configuration du calcul
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER CHOIX
+
+C Affichage du menu de configuration du calcul
+C-----------------------------------------------------------------------
+30 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' MENU DE CONFIGURATION DU CALCUL '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)'Choix 1: parametres de simulation generaux '
+ WRITE(*,*)'Choix 2: parametres du transport solide '
+ WRITE(*,*)'Choix 0: retour au menu principal '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'Votre choix:'
+ READ(*,*)CHOIX
+
+C Ventilation
+C-----------------------------------------------------------------------
+ IF(CHOIX.EQ.1)THEN
+ CALL SAIDONNEE
+ GOTO 30
+ ELSEIF(CHOIX.EQ.2)THEN
+ CALL SAITS
+ GOTO 30
+ ELSEIF(CHOIX.EQ.0)THEN
+ RETURN
+ ENDIF
+ WRITE(*,*) 'MAUVAIS CHOIX, RECOMMENCEZ'
+ GOTO 30
+
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE CONDLI
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER I,NT1,NT2,NT3,NT4,CHOIX,MODFIC
+ DOUBLE PRECISION QYH(100),YH(100),QTH(1000),TH(1000),HIMP(100)
+ & ,TA(100),QSTH(100),D(100),S(100)
+ CHARACTER ETUDE*20,NOMFIC*40
+
+ COMMON/ETUDE/ETUDE
+
+ 30 WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'------------------------------------------------------'
+ WRITE(*,*)'MENU POUR LA SAISIE DES CONDITIONS AUX LIMITES '
+ WRITE(*,*)'------------------------------------------------------'
+ WRITE(*,*)'Choix 1: CL aval Q(Z) ou Z(t) condav.',ETUDE
+ WRITE(*,*)'Choix 2: CL amont Q(t) hydro.',ETUDE
+ WRITE(*,*)'Choix 3: CL amont Z(t) hamont.',ETUDE
+ WRITE(*,*)'Choix 4: CL solide amont Qs(t) hydros.',ETUDE
+ WRITE(*,*)'Choix 0: retour'
+ WRITE(*,*)'------------------------------------------------------'
+ WRITE(*,*)' '
+ WRITE(*,*)'Donnez votre choix:'
+ READ(*,*)CHOIX
+
+C Ventilation
+C-----------------------------------------------------------------------
+ IF(CHOIX.EQ.1)THEN
+ GOTO 201
+ ENDIF
+ IF(CHOIX.EQ.2)THEN
+ GOTO 202
+ ENDIF
+ IF(CHOIX.EQ.3)THEN
+ GOTO 203
+ ENDIF
+ IF(CHOIX.EQ.4)THEN
+ GOTO 204
+ ENDIF
+ IF(CHOIX.EQ.0)THEN
+ RETURN
+ ENDIF
+ WRITE(*,*) 'MAUVAIS CHOIX, RECOMMENCEZ'
+ GOTO 30
+
+C Remplissage du fichier condav.etu (Q(Z) ou Z(t) à l'aval)
+C----------------------------------------------------------
+201 CONTINUE
+ NOMFIC='condav.'//ETUDE
+ OPEN(11,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(*,*)' '
+ WRITE(*,*)'Remplissage du fichier CONDAV.ETUDE contenant'
+ WRITE(*,*)' -> la loi Q(Z) aval si CONDAV=1'
+ WRITE(*,*)' un seul point avec Q=0. si Z constant'
+ WRITE(*,*)' -> la loi Z(t) aval si CONDAV=3'
+ WRITE(*,*)' '
+ WRITE(*,*)' 100 points maximum'
+ WRITE(*,*)' '
+ I=0
+
+101 I=I+1
+ IF(I.EQ.100)THEN
+ WRITE(*,*) 'Dernier couple de valeurs:'
+ ELSE
+ WRITE(*,*) 'Rentrez le ',I,'eme couple (Q,Z) ou (t,Z):'
+ ENDIF
+ READ(*,*) QYH(I),YH(I)
+ WRITE (*,*)'Y a-t-il d''autres valeurs ?'
+ WRITE(*,*) 'Oui--->1 Non--->0'
+ READ(*,*)MODFIC
+ IF(MODFIC.EQ.1)THEN
+ GOTO 101
+ ELSE
+ NT1=I
+ WRITE(11,'(I6)')NT1
+ DO I=1,NT1
+ WRITE(11,'(F15.6,1X,F15.6)')QYH(I),YH(I)
+ ENDDO
+ CLOSE(11)
+ ENDIF
+ GOTO 30
+
+C Remplissage du fichier hydro.etu (loi Q(t) à l'amont)
+C------------------------------------------------------
+202 CONTINUE
+ NOMFIC='hydro.'//ETUDE
+ OPEN(12,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(*,*)' '
+ WRITE(*,*)'Remplissage du fichier HYDRO.ETUDE contenant'
+ WRITE(*,*)' la loi Q(t) amont [CONDAM=1]'
+ WRITE(*,*)' '
+ WRITE(*,*)' 100 points maximum'
+ WRITE(*,*)' '
+ I=0
+
+102 I=I+1
+ IF(I.EQ.1000)THEN
+ WRITE(*,*) 'Dernier couple de valeurs:'
+ ELSE
+ WRITE(*,*) 'Rentrez le ',I,'eme couple (Q,t):'
+ ENDIF
+ READ(*,*) QTH(I),TH(I)
+ WRITE (*,*)'Y a-t-il d''autres valeurs ?'
+ WRITE(*,*) 'Oui--->1 Non--->0'
+ READ(*,*)MODFIC
+ IF(MODFIC.EQ.1)THEN
+ GOTO 102
+ ELSE
+ NT2=I
+ WRITE(12,'(I6)')NT2
+ DO I=1,NT2
+ WRITE(12,'(F15.6,1X,F15.6)')QTH(I),TH(I)
+ ENDDO
+ CLOSE(12)
+ ENDIF
+ GOTO 30
+
+C Remplissage du fichier hamont.etu (loi H(T) à l'amont)
+C-------------------------------------------------------
+203 CONTINUE
+ NOMFIC='hamont.'//ETUDE
+ OPEN(13,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(*,*)' '
+ WRITE(*,*)'Remplissage du fichier HAMONT.ETUDE contenant'
+ WRITE(*,*)' les cotes d''eau amont Z(t) [REGIME=1]'
+ WRITE(*,*)' '
+ WRITE(*,*)' 100 points maximum'
+ WRITE(*,*)' '
+ I=0
+
+103 I=I+1
+ IF(I.EQ.100)THEN
+ WRITE(*,*) 'Dernier couple de valeurs:'
+ ELSE
+ WRITE(*,*) 'Rentrez le ',I,'eme couple (Z,t):'
+ ENDIF
+ READ(*,*) HIMP(I),TA(I)
+ WRITE (*,*)'Y a-t-il d''autres valeurs ?'
+ WRITE(*,*) 'Oui--->1 Non--->0'
+ READ(*,*)MODFIC
+ IF(MODFIC.EQ.1)THEN
+ GOTO 103
+ ELSE
+ NT3=I
+ WRITE(13,'(I6)')NT3
+ DO I=1,NT3
+ WRITE(13,'(F15.6,1X,F15.6)')HIMP(I),TA(I)
+ ENDDO
+ CLOSE(13)
+ ENDIF
+ GOTO 30
+
+C Remplissage du fichier hydros.etu (loi Qs(T) à l'amont)
+C--------------------------------------------------------
+204 CONTINUE
+ NOMFIC='hydros.'//ETUDE
+ OPEN(14,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(*,*)' '
+ WRITE(*,*)'Remplissage du fichier HYDROS.ETUDE contenant'
+ WRITE(*,*)' la loi Qs,D,S(T) amont'
+ WRITE(*,*)' Qs: debit solide (kg/s) si <0, alors capacite solide'
+ WRITE(*,*)' D: diametre reprsentatif (m)'
+ WRITE(*,*)' S: parametre d''etendue (m/m)'
+ WRITE(*,*)' '
+ WRITE(*,*)' 100 points maximum'
+ WRITE(*,*)' '
+ I=0
+
+104 I=I+1
+ IF(I.EQ.100)THEN
+ WRITE(*,*) 'Dernier quadruplet de valeurs:'
+ ELSE
+ WRITE(*,*) 'Rentrez le ',I,'eme quadruplet (Qs,T,D,S):'
+ ENDIF
+ READ(*,*) QSTH(I),TH(I),D(I),S(I)
+ WRITE (*,*)'Y a-t-il d''autres valeurs ?'
+ WRITE(*,*) 'Oui--->1 Non--->0'
+ READ(*,*)MODFIC
+ IF(MODFIC.EQ.1)THEN
+ GOTO 104
+ ELSE
+ NT4=I
+ WRITE(14,'(I3,A,A)')NT4,' Qs(kg/s) temps(s) '
+ & ,'Diam(m) Etend(m/m)'
+ DO I=1,NT4
+ WRITE(14,'(1X,F12.6,1X,G15.6,1X,F12.6,1X,F12.6)')QSTH(I),TH(I)
+ & ,D(I),S(I)
+ ENDDO
+ CLOSE(14)
+ ENDIF
+ GOTO 30
+ END
+
+C-----------------------------------------------------------------------
+ SUBROUTINE DONINI
+C-----------------------------------------------------------------------
+C Saisie des conditions initiales de la simulation
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX,NCMAX,NSMAX
+ PARAMETER(LMAX=3000,NCMAX=300,NSMAX=3000)
+ LOGICAL ENCORE
+ INTEGER CHOIX,I,J,LL,LM,II,N,NBSECT,NUMSEC,IRET,N2,NC
+ DOUBLE PRECISION XSECT(NSMAX)
+ DOUBLE PRECISION XX,XZCOU,XCTDF(LMAX)
+ DOUBLE PRECISION TMAIL(LMAX),XTMAIL(LMAX)
+ DOUBLE PRECISION VINIT(LMAX),VINIT0(NSMAX)
+ DOUBLE PRECISION YINIT(LMAX),YINIT0(NSMAX)
+ DOUBLE PRECISION HINIT(LMAX),HINIT0(NSMAX)
+ DOUBLE PRECISION QSINIT(LMAX),QSINIT0(NSMAX)
+ DOUBLE PRECISION DINIT(LMAX),DINIT0(NSMAX)
+ DOUBLE PRECISION SINIT(LMAX),SINIT0(NSMAX)
+ CHARACTER ETUDE*20,NOMFIC*38,AA*1,REP*1
+ INTEGER IREPVIT
+
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+ COMMON/ETUDE/ETUDE
+
+C Lecture du fichier mail
+C-----------------------------------------------------------------------
+ NOMFIC='mail.'//ETUDE
+ OPEN(20,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ READ(20,*)LM
+ READ(20,*)(TMAIL(I),I=1,LM)
+ READ(20,*)(XTMAIL(I),I=1,LM-1)
+ CLOSE(20)
+
+C Menu
+C-----------------------------------------------------------------------
+ 30 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' MENU DE SAISIE DES CONDITIONS INITIALES '
+ WRITE(*,*)'--------------------------------------------------'
+Cliq WRITE(*,*)'Choix 1: C.I. par zones dans condi0'
+Cliq WRITE(*,*)'Choix 2: interpolation des C.I. en condin'
+ WRITE(*,*)'Choix 1: C.I. liquides par zones dans condi0'
+ WRITE(*,*)'Choix 2: interpolation des C.I. liquides en condin'
+ WRITE(*,*)'Choix 3: C.I. solides par zones dans condis0'
+ WRITE(*,*)'Choix 4: interpolation des C.I. solides en condins'
+ WRITE(*,*)'Choix 0: retour au menu principal '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' '
+ WRITE(*,*)'Votre choix:'
+ READ(*,*)CHOIX
+
+C Ventilation
+C-----------------------------------------------------------------------
+ IF(CHOIX.EQ.1)THEN
+ GOTO 412
+ ELSEIF(CHOIX.EQ.3)THEN
+ GOTO 414
+ ELSEIF(CHOIX.EQ.2)THEN
+ GOTO 413
+ ELSEIF(CHOIX.EQ.4)THEN
+ GOTO 415
+ ELSEIF(CHOIX.EQ.0)THEN
+ RETURN
+ ENDIF
+ WRITE(*,*) 'MAUVAIS CHOIX, RECOMMENCEZ'
+ GOTO 30
+
+C Saisie d'un fichier 'condi0.etude'
+C-----------------------------------------------------------------------
+C Saisie
+412 CONTINUE
+ write(*,*)'voulez vous rentrer '
+ write(*,*)'vitesse ----> 1'
+ write(*,*)'debit ----> 2 '
+ read(*,*)irepvit
+ If(irepvit.NE.2)IREPVIT=1
+ IRET=0
+ ENCORE=.TRUE.
+ NUMSEC=0
+ DO WHILE(ENCORE)
+
+
+ NUMSEC=NUMSEC+1
+ 750 WRITE(*,*)'On est dans la section numero ',NUMSEC
+ WRITE(*,*)'Donnez l''abscisse de la section'
+ WRITE(*,*)'-----> XSECT(',NUMSEC,') = ?'
+ READ (*,*) XSECT(NUMSEC)
+
+C Vérification
+ IF(XSECT(NUMSEC).LT.TMAIL(1).OR.XSECT(NUMSEC).GT.TMAIL(LM))THEN
+ WRITE(*,*)'ATTENTION: Section en dehors du domaine maille'
+ ENDIF
+
+C Entrée des données initiales dans la section
+ WRITE(*,*)'Donnez le tirant d''eau initial (=10000. si retenue)'
+ WRITE(*,*)'YINIT pour l''abscisse',XSECT(NUMSEC),'?'
+ READ(*,*) YINIT0(NUMSEC)
+ IF(YINIT0(NUMSEC).GT.9999.)THEN
+ WRITE(*,*)'Le fichier geom.etude doit etre defini.'
+ WRITE(*,*)'En amont et aval de retenue doivent figurer'
+ WRITE(*,*)'2 sections a la meme abscisse'
+ IF(NUMSEC.NE.1.AND.IRET.EQ.0
+ & .AND.XSECT(NUMSEC).NE.XSECT(NUMSEC-1)) GOTO 750
+ WRITE(*,*)'Donnez la cote initiale dans la retenue'
+ WRITE(*,*)'ZINIT pour l''abscisse',XSECT(NUMSEC)
+ READ(*,*)HINIT0(NUMSEC)
+ IRET=1
+ ELSE
+ IF(IRET.EQ.1.AND.YINIT0(NUMSEC-1).GT.9999.
+ & .AND.XSECT(NUMSEC-1).NE.XSECT(NUMSEC))THEN
+ WRITE(*,*)'ERREUR: abscisse differente de la section precedente'
+ GOTO 750
+ ENDIF
+ ENDIF
+ IF(IREPVIT.EQ.1)THEN
+ WRITE(*,*)'Donnez la vitesse initiale'
+c WRITE(*,*)'Vous pouvez aussi donner le debit initial: dans ce '
+c & ,'cas rentrer ce debit + 1000'
+ WRITE(*,*)'VINIT pour l''abscisse',XSECT(NUMSEC),'?'
+ READ(*,*) VINIT0(NUMSEC)
+ ELSE
+ WRITE(*,*)'Donnez le debit initial'
+ WRITE(*,*)'QINIT pour l''abscisse',XSECT(NUMSEC),'?'
+ READ(*,*) VINIT0(NUMSEC)
+ VINIT0(NUMSEC)=VINIT0(NUMSEC)+1000.
+C fin du if sur vitesse ou debit
+ ENDIF
+ IF(NUMSEC.GT.1)THEN
+ IF(VINIT0(NUMSEC).GT.500..AND.VINIT0(NUMSEC-1).LT.500.)THEN
+ IF(XSECT(NUMSEC).NE.XSECT(NUMSEC-1))THEN
+ WRITE(*,*)'ATTENTION les 2 dernieres abscisses ne sont '
+ & ,'pas egales alors qu''on passe de VINIT a QINIT.'
+ WRITE(*,*)'L''interpolation sera mauvaise'
+ ENDIF
+ ELSEIF(VINIT0(NUMSEC-1).GT.500..AND.VINIT0(NUMSEC).LT.500.)THEN
+ IF(XSECT(NUMSEC).NE.XSECT(NUMSEC-1))THEN
+ WRITE(*,*)'ATTENTION les 2 dernieres abscisses ne sont '
+ & ,'pas egales alors qu''on passe de QINIT a VINIT.'
+ WRITE(*,*)'L''interpolation sera mauvaise'
+ ENDIF
+ ENDIF
+ ENDIF
+
+C Fin de la saisie pour la section
+ WRITE(*,*)'Y a-t-il d''autres sections ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'N').OR.(REP.EQ.'n')) ENCORE=.FALSE.
+ ENDDO
+ NBSECT=NUMSEC
+
+C Ecriture dans 'condi0.etude'
+ WRITE(*,*)' '
+ WRITE(*,*)'Voulez-vous sauvegarder dans un fichier condi0.'
+ & ,ETUDE,' ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+ NOMFIC='condi0.'//ETUDE
+ OPEN(32,FILE=NOMFIC,STATUS='UNKNOWN')
+ DO J=1,NBSECT
+ WRITE(32,*) XSECT(J),YINIT0(J),VINIT0(J),HINIT0(J)
+ ENDDO
+ CLOSE(32)
+ ENDIF
+
+C Interpolation dans 'condin.etude'
+C-----------------------------------------------------------------------
+C Introduction
+ WRITE(*,*)' '
+ WRITE(*,*)'Voulez-vous interpoler dans un fichier condin.'
+ & ,ETUDE,' ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'N').OR.(REP.EQ.'n')) GOTO 30
+413 CONTINUE
+
+C Lecture des données par zone
+ IRET=0
+ NOMFIC='condi0.'//ETUDE
+ OPEN(52,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ NBSECT=0
+ NUMSEC=0
+ DO WHILE(.TRUE.)
+ NUMSEC=NUMSEC+1
+ IF(NUMSEC.GT.NSMAX)THEN
+ WRITE(*,*)'LE FICHIER CONDI0 CONTIENT PLUS DE',NSMAX,'POINTS'
+ WRITE(*,*)'LES POINTS SUPPLEMENTAIRES NE SONT PAS PRIS'
+ GOTO 441
+ ENDIF
+ READ(52,*,END=441) XSECT(NUMSEC),YINIT0(NUMSEC),VINIT0(NUMSEC)
+ & ,HINIT0(NUMSEC)
+ ENDDO
+441 CONTINUE
+ NBSECT=NUMSEC-1
+ CLOSE(52)
+ DO NUMSEC=1,NBSECT
+ IF(YINIT0(NUMSEC).GT.9999.) IRET=1
+ ENDDO
+
+C Interpolation aux centremailles
+ DO J=1,NBSECT-1
+ IF(XSECT(J).GT.XSECT(J+1))THEN
+ WRITE(*,*)'ABSCISSES NON CROISSANTES DANS condi0'
+ STOP
+ ENDIF
+ ENDDO
+ CALL INTERP(XSECT,NBSECT,YINIT,YINIT0)
+ CALL INTERP(XSECT,NBSECT,VINIT,VINIT0)
+ IF(IRET.EQ.1) CALL INTERP(XSECT,NBSECT,HINIT,HINIT0)
+
+C Écriture dans condin.etu
+ N=23
+ NOMFIC='condin.'//ETUDE
+ OPEN(N,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(*,*)LM,' centremailles'
+ WRITE(N,*)'0. ',LM
+ IF(IRET.EQ.1)THEN
+ N2=21
+ NOMFIC='geomac-i.'//ETUDE
+ OPEN(N2,FILE=NOMFIC,STATUS='OLD')
+ READ(N2,*)LL
+ DO J=1,LL
+ READ(N2,*)II,XX,NC
+ XCTDF(J)=9999.
+ DO I=1,NC
+ READ(N2,'(A1,1X,F11.5,1X,F11.5)')AA,XX,XZCOU
+ IF(XZCOU.LT.XCTDF(J)) XCTDF(J)=XZCOU
+ ENDDO
+ ENDDO
+ IF(YINIT(1).GT.9999.)THEN
+ YINIT(1)=HINIT(1)-XCTDF(1)
+ IF(YINIT(1).LT.0.) YINIT(1)=0.
+ ENDIF
+ DO J=2,LM-1
+ IF(YINIT(J).GT.9999.)THEN
+ YINIT(J)=HINIT(J)-.5*(XCTDF(J-1)+XCTDF(J))
+ IF(YINIT(J).LT.0.) YINIT(J)=0.
+ ENDIF
+ ENDDO
+ IF(YINIT(LM).GT.9999.)THEN
+ YINIT(LM)=HINIT(LM)-XCTDF(LL)
+ IF(YINIT(LM).LT.0.) YINIT(LM)=0.
+ ENDIF
+ CLOSE(N2)
+ ENDIF
+ DO I=1,LM
+ WRITE(N,*)I,YINIT(I),VINIT(I)
+ ENDDO
+ CLOSE(N)
+ GOTO 30
+
+C Saisie d'un fichier 'condis0.etude'
+C-----------------------------------------------------------------------
+C Saisie
+414 CONTINUE
+ ENCORE=.TRUE.
+ NUMSEC=0
+ DO WHILE(ENCORE)
+ NUMSEC=NUMSEC+1
+ 751 WRITE(*,*)'On est dans la section numero ',NUMSEC
+ WRITE(*,*)'Donnez l''abscisse de la section'
+ WRITE(*,*)'-----> XSECT(',NUMSEC,') = ?'
+ READ (*,*) XSECT(NUMSEC)
+
+C Vérification
+ IF(XSECT(NUMSEC).LT.TMAIL(1).OR.XSECT(NUMSEC).GT.TMAIL(LM))THEN
+ WRITE(*,*)'ATTENTION: section en dehors du domaine maille'
+ ENDIF
+
+C Entrée des données initiales dans la section
+ WRITE(*,*)'Donnez le debit solide massique (kg/s)'
+ WRITE(*,*)'QSINIT pour l''abscisse',XSECT(NUMSEC),'?'
+ READ(*,*) QSINIT0(NUMSEC)
+c IF(QSINIT0(NUMSEC).LT.0.)THEN
+c WRITE(*,*)'ERREUR: la valeur doit etre positive. Recommencez.
+c GOTO 751
+c ENDIF
+
+ WRITE(*,*)'Donnez le diametre representatif (m)'
+ WRITE(*,*)'DINIT pour l''abscisse',XSECT(NUMSEC),'?'
+ READ(*,*) DINIT0(NUMSEC)
+ IF(DINIT0(NUMSEC).LT.0.)THEN
+ WRITE(*,*)'ERREUR: la valeur doit etre positive. Recommencez.'
+ GOTO 751
+ ENDIF
+
+ WRITE(*,*)'Donnez le parametre d''etendue granulometrique (m/m)'
+ WRITE(*,*)'SINIT pour l''abscisse',XSECT(NUMSEC),'?'
+ READ(*,*) SINIT0(NUMSEC)
+ IF(SINIT0(NUMSEC).LT.0.)THEN
+ WRITE(*,*)'ERREUR: la valeur doit etre positive. Recommencez.'
+ GOTO 751
+ ENDIF
+
+C Fin de la saisie pour la section
+ WRITE(*,*)'Y a-t-il d''autres sections ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'N').OR.(REP.EQ.'n')) ENCORE=.FALSE.
+ ENDDO
+ NBSECT=NUMSEC
+
+C Ecriture dans condis0.etu
+C--------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'Voulez-vous sauvegarder dans un fichier condis0.'
+ & ,ETUDE,' ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+ NOMFIC='condis0.'//ETUDE
+ OPEN(32,FILE=NOMFIC,STATUS='UNKNOWN')
+ DO J=1,NBSECT
+ WRITE(32,*)XSECT(J),QSINIT0(J),DINIT0(J),SINIT0(J)
+ & ,' X(m) Qs(kg/s) D(m) S(m/m)'
+ ENDDO
+ CLOSE(32)
+ ENDIF
+
+C Interpolation dans 'condins.etude'
+C-----------------------------------------------------------------------
+C Introduction
+ WRITE(*,*)' '
+ WRITE(*,*)'Voulez-vous interpoler dans un fichier condins.'
+ & ,ETUDE,' ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'N').OR.(REP.EQ.'n')) GOTO 30
+415 CONTINUE
+
+C Lecture des données par zone
+ N=52
+ NOMFIC='condis0.'//ETUDE
+ OPEN(N,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ NBSECT=0
+ NUMSEC=0
+ DO WHILE(.TRUE.)
+ NUMSEC=NUMSEC+1
+ IF(NUMSEC.GT.NSMAX)THEN
+ WRITE(*,*)'LE FICHIER CONDIS0 CONTIENT PLUS DE',NSMAX,'POINTS'
+ WRITE(*,*)'LES POINTS SUPPLEMENTAIRES NE SONT PAS PRIS'
+ GOTO 443
+ ENDIF
+ READ(N,*,END=443)XSECT(NUMSEC),QSINIT0(NUMSEC),DINIT0(NUMSEC)
+ & ,SINIT0(NUMSEC)
+ ENDDO
+443 NBSECT=NUMSEC-1
+ CLOSE(N)
+
+C Interpolation aux intermailles
+ DO J=1,NBSECT-1
+ IF(XSECT(J).GT.XSECT(J+1))THEN
+ WRITE(*,*)'ABSCISSES NON CROISSANTES DANS condis0'
+ STOP
+ ENDIF
+ ENDDO
+ CALL INTER2(XSECT,NBSECT,QSINIT,QSINIT0)
+ CALL INTER2(XSECT,NBSECT,DINIT,DINIT0)
+ CALL INTER2(XSECT,NBSECT,SINIT,SINIT0)
+
+C Ecriture dans condins.etu
+ N=23
+ NOMFIC='condins.'//ETUDE
+ OPEN(N,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(*,*)LM-1,' intermailles'
+ WRITE (N,*)' 0. NoIM Qs(kg/s) D(m) S(m/m)'
+ DO I=1,LM-1
+ WRITE(N,*)I,QSINIT(I),DINIT(I),SINIT(I)
+ ENDDO
+ CLOSE(N)
+ GOTO 30
+
+C Traitement des erreurs
+C-----------------------------------------------------------------------
+ 1000 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'ERREUR sur le fichier ',NOMFIC
+ WRITE(*,*)'Le fichier n''existe pas ou est illisible'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+ RETURN
+ END
+
+
+
+C-----------------------------------------------------------C
+ SUBROUTINE DXCONS
+C-----------------------------------------------------------C
+C SON ROLE : selection d'un maillage de pas constant C
+C-----------------------------------------------------------C
+C IMPLICIT NONE
+ INTEGER LMAX
+ PARAMETER (LMAX=3000)
+
+ INTEGER I,LM,CHOIXM,N
+ CHARACTER*38 NOMFIC,ETUDE*20
+ DOUBLE PRECISION DX,TMAIL(LMAX),XZERO,XTMAIL(LMAX)
+
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+ COMMON/ETUDE/ETUDE
+
+ 5000 WRITE (*,*) 'DONNEZ LE PAS D ESPACE DX'
+ WRITE (*,*) '----> DX = ?'
+ READ (*,*) DX
+
+ WRITE (*,*) 'DONNEZ LE NOMBRE DE POINTS DE CALCUL LM'
+ WRITE (*,*) '( LA DISTANCE TOTALE EST:(LM-2)*DX )'
+ WRITE (*,*) '----> LM = ?'
+ READ (*,*) LM
+ IF(LM.GT.LMAX)THEN
+ WRITE(*,*)'LE NOMBRE DE MAILLES EST LIMITE A',LMAX
+ GO TO 5000
+ ENDIF
+
+C -- SAISIE DE L'ABSCISSE INITIALE
+ WRITE (*,*) 'DONNEZ L ABSCISSE DU PREMIER POINT'
+ WRITE (*,*) '----> XZERO = ?'
+ READ (*,*) XZERO
+
+C -- DOUBLE PRECISIONISATION DU TABLEAU TMAIL DE DIMENSION LM ET CONTENA
+C ---- LES ABSCISSES DES POINTS DE CALCUL
+
+ TMAIL(1)=XZERO
+ XTMAIL(1)=XZERO
+ DO 10 I=2,LM-1
+ TMAIL(I)=DX*(FLOAT(I)-1.5)+XZERO
+ XTMAIL(I)=DX*(FLOAT(I)-1.)+XZERO
+10 CONTINUE
+ TMAIL(LM)=XTMAIL(LM-1)
+C TMAIL(LM)=FLOAT(LM-2)*DX+XZERO
+ WRITE(*,*)'LE DERNIER POINT EST D''ABSCISSE',TMAIL(LM)
+
+C -- ECRITURE DU FICHIER MAIL0.ETUDE
+C---------------------------------------
+
+ WRITE(*,*) 'VOULEZ VOUS SAUVEGARDER DANS UN FICHIER mail0.'
+ +,ETUDE,'?'
+ WRITE(*,*) ' OUI --> 1 NON --> 0'
+ READ (*,*) CHOIXM
+
+ IF (CHOIXM.EQ.1) THEN
+ N=30
+ NOMFIC='mail0.'//ETUDE
+ OPEN(N,FILE=NOMFIC,STATUS='UNKNOWN')
+
+ WRITE(N,*) XZERO,DX,DX
+ WRITE(N,*) TMAIL(LM),DX,DX
+
+ CLOSE(N)
+ ENDIF
+
+C -- ECRITURE DU FICHIER MAIL.ETUDE
+ CALL ECRIMA
+
+ END
+
+
+C-------------------------------------------------------------C
+ SUBROUTINE DXVARE
+C-------------------------------------------------------------C
+C SON ROLE : selection d'un maillage à pas variables C
+C de pas d'espace variant lineairement entre C
+C deux sections, une de debut et une de fin C
+C-------------------------------------------------------------C
+C IMPLICIT NONE
+ INTEGER LMAX,NSMAX
+
+ PARAMETER (LMAX=3000,NSMAX=3000)
+
+ INTEGER ENCORE,N,NBSECT,NPAS,NUMSEC
+
+ DOUBLE PRECISION DXSECT(NSMAX),LONGSE,DXSEC2(NSMAX)
+ DOUBLE PRECISION XSECT(NSMAX)
+
+ CHARACTER ETUDE*20,NOMFIC*38
+
+ COMMON/ETUDE/ETUDE
+
+ N=30
+ NOMFIC='mail0.'//ETUDE
+ OPEN(N,FILE=NOMFIC,STATUS='UNKNOWN')
+C 1)-->ENTREE DES SECTIONS ET DU PAS D'ESPACE CORRESPONDANT
+C---------------------------------------------------------------
+C -------------- NBSECT = NOMBRE DE SECTEURS ------------------C
+C--------------- NUMSEC = NUMERO DU SECTEUR -------------------C
+
+ WRITE (*,*) 'DONNEZ L ABSCISSE DE DEPART'
+ WRITE (*,*) '----> XSECT(1) = ?'
+ READ (*,*) XSECT(1)
+ NUMSEC=0
+10 NUMSEC=NUMSEC+1
+ WRITE (*,*) 'ON EST DANS LE SECTEUR NUMERO ',NUMSEC
+ WRITE (*,*) 'DONNEZ L ABSCISSE DE L EXTREMITE AVAL DU SECTEUR'
+ WRITE (*,*) '----> XSECT(',NUMSEC + 1,') = ?'
+ READ (*,*) XSECT(NUMSEC+1)
+
+C ---- VERIFICATION ABSCISSE
+ IF(XSECT(NUMSEC+1).LE.XSECT(NUMSEC)) THEN
+ WRITE(*,*) 'XSECT EST TROP PETIT'
+ WRITE(*,*) 'RECOMMENCEZ'
+ NUMSEC=NUMSEC-1
+ GOTO 10
+ ENDIF
+
+20 WRITE(*,*) 'DONNEZ LE PAS D ESPACE DX DESIRE AMONT SECTEUR'
+ WRITE(*,*) '----> DXSECT(',NUMSEC,') = ?'
+ READ (*,*) DXSECT(NUMSEC)
+C ---- VERIFICATION PAS D'ESPACE
+ IF(DXSECT(NUMSEC).GT.(XSECT(NUMSEC+1)-XSECT(NUMSEC))) THEN
+ WRITE(*,*) 'DXSECT EST TROP GRAND'
+ WRITE (*,*) 'RECOMMENCEZ'
+ GOTO 20
+ ENDIF
+ 21 WRITE(*,*) 'DONNEZ LE PAS D ESPACE DX DESIRE AVAL SECTEUR'
+ WRITE(*,*) '----> DXSEC2(',NUMSEC,') = ?'
+ READ (*,*) DXSEC2(NUMSEC)
+C ---- VERIFICATION PAS D'ESPACE
+ IF(DXSEC2(NUMSEC).GT.(XSECT(NUMSEC+1)-XSECT(NUMSEC))) THEN
+ WRITE(*,*) 'DXSEC2 EST TROP GRAND'
+ WRITE (*,*) 'RECOMMENCEZ'
+ GOTO 21
+ ENDIF
+ WRITE(N,*) XSECT(NUMSEC),DXSECT(NUMSEC),DXSEC2(NUMSEC)
+
+
+11 WRITE(*,*) 'Y A T IL D AUTRES SECTEURS ?'
+ WRITE(*,*) ' OUI --> 1 NON --> 0'
+ READ (*,*) ENCORE
+ IF (ENCORE.EQ.1)THEN
+ GOTO 10
+ ELSE IF(ENCORE.NE.0)THEN
+ WRITE(*,*) 'ERREUR DE REPONSE'
+ GO TO 11
+ ENDIF
+ NBSECT=NUMSEC
+ WRITE(N,'(3F15.3)')XSECT(NBSECT+1),DXSECT(NBSECT),DXSEC2(NBSECT)
+
+C 2)--> ON CALCULE LE NOMBRE DE MAILLES ENTIER
+C---------------------------------------------------------------
+
+ WRITE(*,*) 'VOUS AVEZ RENTRE ',NBSECT,' SECTEURS'
+ WRITE(*,*) 'CALCUL DU NOMBRE DE MAILLES'
+
+ DO 30 NUMSEC=1,NBSECT
+ LONGSE=XSECT(NUMSEC+1)-XSECT(NUMSEC)
+ NPAS=NINT(2.*LONGSE/(DXSECT(NUMSEC)+DXSEC2(NUMSEC)))
+C DXSECT(NUMSEC)=LONGSE/FLOAT(NPAS)
+ WRITE(*,*) 'SECTEUR NUMERO',NUMSEC,
+ :' DEBUT EN',XSECT(NUMSEC),' FIN EN',XSECT(NUMSEC+1),
+ : ' PAS D''ESPACE MOYEN ',LONGSE/FLOAT(NPAS),
+ : ' NOMBRE DE MAILLES ',NPAS
+30 CONTINUE
+
+C 3)--> ECRITURE DU FICHIER MAIL0.ETUDE
+C---------------------------------------
+
+C WRITE(*,*) 'VOULEZ VOUS SAUVEGARDER DANS UN FICHIER MAIL0.'
+C +,ETUDE,'?'
+C WRITE(*,*) ' OUI --> 1 NON --> 0'
+C READ (*,*) CHOIXM
+
+
+C IF (CHOIXM.EQ.1) THEN
+C N=30
+C NOMFIC='MAIL0.'//ETUDE
+C OPEN(N,FILE=NOMFIC,STATUS='UNKNOWN')
+C DO 40 J=1,NBSECT
+C WRITE(N,*) XSECT(J),DXSECT(J),DXSEC2(J)
+C40 CONTINUE
+C WRITE(N,*) XSECT(NBSECT+1),DXSECT(NBSECT),DXSEC2(NBSECT)
+
+ CLOSE(N)
+C ENDIF
+
+C 4)--> DOUBLE PRECISIONISATION DU MAILLAGE
+C--------------------------------
+
+ CALL RELTMA(XSECT,DXSECT,DXSEC2,NBSECT)
+ CALL ECRIMA
+
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE DXVARL
+C-----------------------------------------------------------------------
+C SON ROLE : le fichier de maillage existant il ne reste plus
+C qu'à interpoler
+C-----------------------------------------------------------------------
+C IMPLICIT NONE
+ INTEGER LMAX,NSMAX
+
+ PARAMETER (LMAX=3000,NSMAX=3000)
+
+ INTEGER N,NBSECT,NUMSEC
+
+ DOUBLE PRECISION DXSECT(NSMAX),DXSEC2(NSMAX)
+ DOUBLE PRECISION XSECT(NSMAX)
+
+ CHARACTER ETUDE*20,NOMFIC*38
+
+ COMMON/ETUDE/ETUDE
+
+C 1)--> LECTURE DE MAIL0.ETUDE ET RECUPERATION DES VARIABLES
+C------------------------------------------------------------
+
+ N=30
+ NOMFIC='mail0.'//ETUDE
+ OPEN(N,FILE=NOMFIC,STATUS='OLD')
+
+ NBSECT=0
+ NUMSEC=0
+10 NUMSEC=NUMSEC+1
+C READ(N,'(A80)',END=20)CHAIN
+C READ(CHAIN,'(3F15.3)')XSECT(NUMSEC),DXSECT(NUMSEC),DXSEC2(NUMSEC
+ READ(N,*,END=20)XSECT(NUMSEC),DXSECT(NUMSEC),DXSEC2(NUMSEC)
+C IF(DXSEC2(NUMSEC).EQ.0.)DXSEC2(NUMSEC)=DXSECT(NUMSEC)
+ GOTO 10
+20 NBSECT=NUMSEC-2
+
+ CLOSE(N)
+
+C 2)-->INTERPOLATION
+C---------------------
+
+ CALL RELTMA(XSECT,DXSECT,DXSEC2,NBSECT)
+ CALL ECRIMA
+
+ END
+
+
+C------------------------------------------------------------------C
+ SUBROUTINE ECRIMA
+C------------------------------------------------------------------C
+C SON ROLE : ecriture du maillage dans le fichier mail.etude C
+C------------------------------------------------------------------C
+C IMPLICIT NONE
+ INTEGER LMAX
+ PARAMETER(LMAX=3000)
+
+ INTEGER I,LM,N
+ DOUBLE PRECISION TMAIL(LMAX),XTMAIL(LMAX)
+ CHARACTER ETUDE*20,NOMFIC*38
+
+ COMMON/ETUDE/ETUDE
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+10 N=20
+ IF(LM.GT.LMAX)THEN
+ WRITE(*,*)'LE NOMBRE DE MAILLES EST LIMITE A',LMAX
+ WRITE(*,*)'RECOMMENCEZ LA SAISIE DU MAILLAGE'
+ RETURN
+ ENDIF
+ NOMFIC='mail.'//ETUDE
+ OPEN(N,FILE=NOMFIC,STATUS='UNKNOWN')
+C WRITE (*,*) LM
+C WRITE (*,*) (TMAIL(I),I=1,LM)
+ WRITE (N,*) LM
+ WRITE (N,*) (TMAIL(I),I=1,LM)
+ WRITE (N,*) (XTMAIL(I),I=1,LM-1)
+ CLOSE(N)
+ RETURN
+ END
+
+C-----------------------------------------------------------------------
+ SUBROUTINE FROTTE
+C-----------------------------------------------------------------------
+C Menu relatif aux frottements
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER CHOIX
+ CHARACTER ETUDE*20,NOMFIC0*40,NOMFIC*40,TYPLIT*10
+
+ COMMON/ETUDE/ETUDE
+
+C Affichage du menu
+C-----------------------------------------------------------------------
+30 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' MENU DE SAISIE DES FROTTEMENTS (STRICKLER) '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)'Choix 1: saisie par zones des coefficients de '
+ WRITE(*,*)' Strickler du lit mineur dans frot0 '
+ WRITE(*,*)'Choix 2: saisie par zones des coefficients de '
+ WRITE(*,*)' Strickler du lit moyen dans frot20 '
+ WRITE(*,*)'Choix 3: interpolation des coefficients de '
+ WRITE(*,*)' Strickler du lit mineur dans frot '
+ WRITE(*,*)'Choix 4: interpolation des coefficients de '
+ WRITE(*,*)' Strickler du lit moyen dans frot2 '
+ WRITE(*,*)'Choix 0: retour au menu principal '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'Votre choix:'
+ READ(*,*)CHOIX
+
+C Ventilation
+C-----------------------------------------------------------------------
+ IF(CHOIX.EQ.1)THEN
+ NOMFIC0='frot0.'//ETUDE
+ NOMFIC='frot.'//ETUDE
+ TYPLIT='lit mineur'
+ CALL FROTT0(TYPLIT,NOMFIC0,NOMFIC)
+ GOTO 30
+ ELSEIF(CHOIX.EQ.2)THEN
+ NOMFIC0='frot20.'//ETUDE
+ NOMFIC='frot2.'//ETUDE
+ TYPLIT='lit moyen'
+ CALL FROTT0(TYPLIT,NOMFIC0,NOMFIC)
+ GOTO 30
+ ELSEIF(CHOIX.EQ.3)THEN
+ NOMFIC0='frot0.'//ETUDE
+ NOMFIC='frot.'//ETUDE
+ TYPLIT='lit mineur'
+ CALL FROTT(TYPLIT,NOMFIC0,NOMFIC)
+ GOTO 30
+ ELSEIF(CHOIX.EQ.4)THEN
+ NOMFIC0='frot20.'//ETUDE
+ NOMFIC='frot2.'//ETUDE
+ TYPLIT='lit moyen'
+ CALL FROTT(TYPLIT,NOMFIC0,NOMFIC)
+ GOTO 30
+ ELSEIF(CHOIX.EQ.0)THEN
+ RETURN
+ ENDIF
+ WRITE(*,*) 'MAUVAIS CHOIX, RECOMMENCEZ'
+ GOTO 30
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE FROTT0(TYPLIT,NOMFIC0,NOMFIC)
+C-----------------------------------------------------------------------
+C Saisie des frottements par zones
+C-----------------------------------------------------------------------
+C TYPLIT indique le type de lit (mineur ou moyen)
+C NOMFIC0 le nom du fichier de sauvegarde de la saisie par zones
+C NOMFIC le nom du fichier de sauvegarde des données interpolées
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER NSMAX
+ PARAMETER (NSMAX=3000)
+ LOGICAL ENCORE
+ INTEGER I,NBSECT
+ DOUBLE PRECISION XSECT(NSMAX),FROT0(NSMAX)
+ CHARACTER TYPLIT*10,NOMFIC0*40,NOMFIC*40,REP*1
+
+C Initialisation
+C-----------------------------------------------------------------------
+ ENCORE=.TRUE.
+ I=0
+ WRITE(*,*)' '
+ WRITE(*,*)'SAISIE PAR ZONES DES COEFFICIENTS DE STRICKLER'
+ WRITE(*,*)'On est dans le ',TYPLIT
+ WRITE(*,*)' '
+
+C Saisie des largeurs actives par zones
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'ATTENTION: Les abscisses doivent etre croissantes.'
+ DO WHILE(ENCORE)
+ I=I+1
+ WRITE(*,*)' '
+ 1001 CONTINUE
+ WRITE(*,*)'Entrez le couple abscisse (m) - Strickler (m(1/3)/s)'
+ & ,' numero ',I
+ READ(*,*)XSECT(I),FROT0(I)
+ IF((I.GT.1).AND.(XSECT(I).LE.XSECT(I-1)))THEN
+ WRITE(*,*)'ATTENTION: abscisse non croissante. Recommencez.'
+ GOTO 1001
+ ENDIF
+ IF(FROT0(I).LT.0.)THEN
+ WRITE(*,*)'ATTENTION: le Strickler ne peut pas etre negatif. '
+ & ,'Recommencez.'
+ GOTO 1001
+ ENDIF
+ WRITE(*,*)'Y-a-il un autre couple ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'N').OR.(REP.EQ.'n')) ENCORE=.FALSE.
+ ENDDO
+ NBSECT=I
+
+C Écriture du fichier de données par zones
+C-----------------------------------------------------------------------
+ OPEN(20,FILE=NOMFIC0,STATUS='UNKNOWN')
+ DO I=1,NBSECT
+ WRITE(20,120)XSECT(I),FROT0(I)
+ ENDDO
+ CLOSE(20)
+ WRITE(*,*)' '
+ WRITE(*,*)'Ecriture du fichier ',NOMFIC0
+
+C Interpolation éventuelle
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'Voulez-vous interpoler pour le ',TYPLIT,' ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o')) CALL FROTT(TYPLIT,NOMFIC0,NOMFIC)
+
+ RETURN
+ 120 FORMAT(1X,F11.3,1X,F11.5)
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE FROTT(TYPLIT,NOMFIC0,NOMFIC1)
+C-----------------------------------------------------------------------
+C Interpolation des frottements
+C-----------------------------------------------------------------------
+C TYPLIT indique le type de lit (mineur ou moyen)
+C NOMFIC0 le nom du fichier de sauvegarde de la saisie par zones
+C NOMFIC le nom du fichier de sauvegarde des données interpolées
+ IMPLICIT NONE
+ INTEGER LMAX,NSMAX
+ PARAMETER(LMAX=3000,NSMAX=3000)
+ INTEGER I,LM,NBSECT,NUMSEC
+ DOUBLE PRECISION XSECT(NSMAX),FROT(LMAX),FROT0(NSMAX)
+ DOUBLE PRECISION TMAIL(LMAX),XTMAIL(LMAX)
+ CHARACTER ETUDE*20,TYPLIT*10,NOMFIC*40,NOMFIC0*40,NOMFIC1*40
+
+ COMMON/ETUDE/ETUDE
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+C Lecture du fichier de maillage
+C-----------------------------------------------------------------------
+ NOMFIC='mail.'//ETUDE
+ OPEN(20,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ READ(20,*)LM
+ READ(20,*)(TMAIL(I),I=1,LM)
+ READ(20,*)(XTMAIL(I),I=1,LM-1)
+ CLOSE(20)
+
+C Lecture du fichier de frottements par zones
+C-----------------------------------------------------------------------
+ NOMFIC=NOMFIC0
+ OPEN(30,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ NUMSEC=1
+ DO WHILE(.TRUE.)
+ READ(30,*,END=300)XSECT(NUMSEC),FROT0(NUMSEC)
+ NUMSEC=NUMSEC+1
+ ENDDO
+ 300 CONTINUE
+ NBSECT=NUMSEC-1
+ CLOSE(30)
+
+C Interpolation
+C-----------------------------------------------------------------------
+ CALL INTERP(XSECT,NBSECT,FROT,FROT0)
+
+C Écriture du fichier interpolé
+C-----------------------------------------------------------------------
+ NOMFIC=NOMFIC1
+ OPEN(20,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(20,131)LM
+ DO I=1,LM
+ WRITE(20,132)I,FROT(I)
+ ENDDO
+ CLOSE(20)
+ WRITE(*,*)' '
+ WRITE(*,*)'Pour le ',TYPLIT,', ecriture du fichier ',NOMFIC
+ RETURN
+
+C Traitement des erreurs
+C-----------------------------------------------------------------------
+ 1000 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'ERREUR sur le fichier ',NOMFIC
+ WRITE(*,*)'Le fichier n''existe pas ou est illisible'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+ RETURN
+
+ 131 FORMAT(1X,I4)
+ 132 FORMAT(1X,I4,1X,F11.5)
+ END
+
+
+C------------------------------------------------------------------
+ SUBROUTINE INTERP(XS,NS,Y,Y0)
+C------------------------------------------------------------------
+C SON ROLE : interpoler y0 (-->y) entre deux sections C
+C correspondant aux points de maillage XS C
+C entrees : C
+C C
+C ---- LM EST LE NOMBRE DE MAILLE C
+C ---- NS LE NOMBRE DES SECTIONS DE DONNEE C
+C ---- TMAIL LES ABSCISSES DES POINTS DU MAILLAGE C
+C ---- XS LES ABSCISSES DES SECTIONS DE DONNEE C
+C ---- Y0 VALEUR AUX SECTIONS DE LA VARIABLE A INTERPOLER C
+C C
+C sorties : C
+C C
+C ---- Y VALEUR AUX POINTS DE MAILLAGE DE LA VARIABLE C
+C A INTERPOLER C
+C C
+C-----------------------------------------------------------------C
+C IMPLICIT NONE
+
+ INTEGER LMAX,NSMAX
+ PARAMETER(LMAX=3000,NSMAX=3000)
+
+ INTEGER J,LM,N,NS,N0,N3
+
+ DOUBLE PRECISION TMAIL(LMAX),XS(NSMAX),Y(LMAX),Y0(NSMAX)
+ :,XTMAIL(LMAX)
+
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+ IF (NS.EQ.1)THEN
+ DO 7 J=1,LM
+ Y(J)=Y0(1)
+7 CONTINUE
+ ELSE
+ N3=1
+ DO 8 J=1,LM
+ N0=N3
+ DO 9 N=N0,NS
+ IF (TMAIL(J).LT.XS(N)) THEN
+ IF(N.EQ.1)THEN
+ Y(J)=Y0(1)
+ ELSE
+ Y(J)=Y0(N-1)+(TMAIL(J)-XS(N-1))*(Y0(N)-Y0(N-1))/(XS(N)-XS(N-1))
+ ENDIF
+ N3=N
+ GO TO 8
+ ENDIF
+9 CONTINUE
+ Y(J)=Y0(NS)
+ N3=NS+1
+8 CONTINUE
+
+
+ ENDIF
+ RETURN
+ END
+
+
+
+C------------------------------------------------------------------
+ SUBROUTINE INTER2(XS,NS,Y,Y0)
+C------------------------------------------------------------------
+C SON ROLE : interpoler y0 (-->y) entre deux sections C
+C correspondant aux points de maillage XS C
+C entrees : C
+C C
+C ---- LM EST LE NOMBRE DE MAILLES C
+C ---- NS LE NOMBRE DES SECTIONS DE DONNEE C
+C ---- XTMAIL LES ABSCISSES DES INTERMAILLES C
+C ---- XS LES ABSCISSES DES SECTIONS DE DONNEE C
+C ---- Y0 VALEUR AUX SECTIONS DE LA VARIABLE A INTERPOLER C
+C C
+C sorties : C
+C C
+C ---- Y VALEUR DE LA VARIABLE A INTERPOLER C
+C AUX POINTS DE MAILLAGE (INTERMAILLES) C
+C C
+C-----------------------------------------------------------------C
+C IMPLICIT NONE
+
+ INTEGER LMAX,NSMAX
+ PARAMETER(LMAX=3000,NSMAX=3000)
+
+ INTEGER J,LM,N,NS,N0,N3
+
+ DOUBLE PRECISION TMAIL(LMAX),XS(NSMAX),Y(LMAX),Y0(NSMAX)
+ :,XTMAIL(LMAX)
+
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+ IF (NS.EQ.1)THEN
+ DO 7 J=1,LM-1
+ Y(J)=Y0(1)
+7 CONTINUE
+ ELSE
+ N3=1
+ DO 8 J=1,LM-1
+ N0=N3
+ DO 9 N=N0,NS
+ IF (XTMAIL(J).LT.XS(N)) THEN
+ IF(N.EQ.1)THEN
+ Y(J)=Y0(1)
+ ELSE
+ Y(J)=Y0(N-1)+(XTMAIL(J)-XS(N-1))*(Y0(N)-Y0(N-1))/(XS(N)-XS(N-1))
+ ENDIF
+ N3=N
+ GO TO 8
+ ENDIF
+9 CONTINUE
+ Y(J)=Y0(NS)
+ N3=NS+1
+8 CONTINUE
+
+
+ ENDIF
+ RETURN
+ END
+
+
+
+C------------------------------------------------------------------
+ SUBROUTINE INTER3(XS,NS,Y,Y0)
+C------------------------------------------------------------------
+C SON ROLE : interpoler y0 (-->y) entre deux sections C
+C correspondant aux points de maillage XS C
+C INTERPOLATION LISSEE (SPLINE) AU LIEU DE LINEAIRE
+C entrees : C
+C C
+C ---- LM EST LE NOMBRE DE MAILLES C
+C ---- NS LE NOMBRE DES SECTIONS DE DONNEE C
+C ---- XTMAIL LES ABSCISSES DES INTERMAILLES C
+C ---- XS LES ABSCISSES DES SECTIONS DE DONNEE C
+C ---- Y0 VALEUR AUX SECTIONS DE LA VARIABLE A INTERPOLER C
+C C
+C sorties : C
+C C
+C ---- Y VALEUR DE LA VARIABLE A INTERPOLER C
+C AUX POINTS DE MAILLAGE (INTERMAILLES) C
+C C
+C-----------------------------------------------------------------C
+C IMPLICIT NONE
+
+ INTEGER LMAX,NSMAX
+ PARAMETER(LMAX=3000,NSMAX=3000)
+
+ INTEGER J,LM,N,NS,N0,N3
+
+ DOUBLE PRECISION TMAIL(LMAX),XS(NSMAX),Y(LMAX),Y0(NSMAX)
+ :,XTMAIL(LMAX)
+ DOUBLE PRECISION H(LMAX),U(LMAX),A,B,C,H1(LMAX),H2(LMAX)
+ :,H3(LMAX)
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+ IF (NS.EQ.1)THEN
+ DO 7 J=1,LM-1
+ Y(J)=Y0(1)
+7 CONTINUE
+ ELSE
+C CALCUL DES COEFFICIENTS DES FONCTIONS SPLINES
+ DO 11 N=1,NS-1
+ H(N)=XS(N+1)-XS(N)
+ 11 CONTINUE
+ DO 12 N=1,NS-2
+ H1(N)=(H(N)+H(N+1))/3.
+ H2(N)=H(N+1)/6.
+ H3(N)=(Y0(N+2)-Y0(N+1))/H(N+1)-(Y0(N+1)-Y0(N))/H(N)
+ 12 CONTINUE
+ DO 13 N=2,NS-2
+ C=H2(N-1)/H1(N-1)
+ H1(N)=H1(N)-H2(N-1)*C
+ H3(N)=H3(N)-H3(N-1)*C
+ 13 CONTINUE
+ U(NS-1)=H3(NS-2)/H1(NS-2)
+ U(1)=0.
+ U(NS)=0.
+ DO 14 N=NS-2,2,-1
+ U(N)=(H3(N-1)-H2(N-1)*U(N+1))/H1(N-1)
+ 14 CONTINUE
+C UTILISATION DES FONCTIONS SPLINES
+ N3=1
+ DO 8 J=1,LM-1
+ N0=N3
+ DO 9 N=N0,NS
+ IF (XTMAIL(J).LT.XS(N)) THEN
+ IF(N.EQ.1)THEN
+ Y(J)=Y0(1)
+ ELSE
+C Y(J)=Y0(N-1)+(XTMAIL(J)-XS(N-1))*(Y0(N)-Y0(N-1))/(XS(N)-XS(N-1))
+ A=(XS(N)-XTMAIL(J))/H(N-1)
+ B=1.-A
+ Y(J)=H(N-1)**2*(A*(A**2-1.)*U(N-1)+B*(B**2-1.)*U(N))/6.
+ Y(J)=Y(J)+A*Y0(N-1)+B*Y0(N)
+ ENDIF
+ N3=N
+ GO TO 8
+ ENDIF
+9 CONTINUE
+ Y(J)=Y0(NS)
+ N3=NS+1
+8 CONTINUE
+
+
+ ENDIF
+ RETURN
+ END
+
+
+
+C-----------------------------------------------------------------------
+
+
+C----------------------------------------------------------------C
+ SUBROUTINE MAILLA
+C----------------------------------------------------------------C
+C SON ROLE : saisie du maillage C
+C----------------------------------------------------------------C
+C IMPLICIT NONE
+
+C PARAMETER (LMAX=3000,NSMAX=3000)
+
+ INTEGER CHOIXM,I
+C DOUBLE PRECISION TMAIL(LMAX),XTMAIL(LMAX)
+
+ CHARACTER ETUDE*20
+ LOGICAL TESTCH
+
+ COMMON/ETUDE/ETUDE
+
+11 WRITE(*,*) ' '
+ WRITE(*,*) ' '
+ WRITE(*,*) '---------------------------------------------------'
+ WRITE(*,*) ' MENU D APPEL DES MAILLEURS'
+ WRITE(*,*) '---------------------------------------------------'
+ WRITE(*,*) 'CHOIX 1 : MAILLAGE DE PAS CONSTANT'
+ WRITE(*,*) 'CHOIX 2 : MAILLAGE VARIABLE ET ECRITURE D UN',
+ :' FICHIER mail0.', ETUDE
+ WRITE(*,*) 'CHOIX 3 : MAILLAGE VARIABLE PAR LECTURE D UN',
+ :' FICHIER mail0.',ETUDE
+ WRITE(*,*) '---------------------------------------------------'
+ WRITE(*,*) ' '
+ WRITE(*,*) 'DONNEZ VOTRE CHOIX'
+ READ (*,*) CHOIXM
+
+C ---- TEST DU CHOIX
+
+ TESTCH=.FALSE.
+ DO 101 I=1,3
+ TESTCH=TESTCH.OR.(CHOIXM.EQ.I)
+101 CONTINUE
+ IF (.NOT.TESTCH) THEN
+ WRITE(*,*) 'VOUS AVEZ DONNE UN MAUVAIS CHOIX'
+ WRITE(*,*) 'RECOMMENCEZ'
+ GOTO 11
+ ENDIF
+
+ GOTO (12,13,14),CHOIXM
+
+C 1)--> MAILLEUR A PAS CONSTANT
+12 CALL DXCONS
+ RETURN
+
+C 2)--> MAILLEUR PAS VARIABLE ET ECRITURE FICHIER MAIL0.ETUDE
+13 CALL DXVARE
+ RETURN
+
+C 3)--> MAILLEUR PAS VARIABLE PAR LECTURE FICHIER MAIL0.ETUDE
+14 CALL DXVARL
+ RETURN
+
+ END
+
+
+C------------------------------------------------------------------C
+ SUBROUTINE RELTMA(XSECT,DXSECT,DXSEC2,NBSECT)
+C------------------------------------------------------------------C
+C SON ROLE : DOUBLE PRECISIONiser un tableau TMAIL de LM elements
+C contenant les abscisses des points a C
+C calculer C
+C C
+C ---- ON DEFINIT LES CARACTERISTIQUES DU MAILLAGE DANS LE SECTEUR C
+C NUMSEC C
+C ---- DXE = TAILLE EFFECTIVE DE LA MAILLE C
+C ---- DXS = TAILLE DE MAILLE SOUHAITEE DANS NUMSEC C
+C ---- JDEBUT = INDICE DE LA PREMIERE INTERFACE A CALCULER C
+C ---- JFIN = INDICE DU DERNIER MILIEU DE MAILLE CALCULE C
+C ---- LMSEC = NOMBRE DE MAILLES DANS LE SECTEUR C
+C ---- XDEBUT = ABSCISSE DU DEBUT C
+C ---- XFIN = ABSCISSE DE FIN C
+C------------------------------------------------------------------C
+C IMPLICIT NONE
+ INTEGER LMAX,NSMAX
+
+ PARAMETER (LMAX=3000,NSMAX=3000)
+
+ INTEGER I,J,JDEBUT,JFIN,LM,NBSECT,NUMSEC,LMSEC
+
+ DOUBLE PRECISION DXS,DXSECT(NSMAX),TMAIL(LMAX),DXE,DXSEC2(NSMAX)
+ DOUBLE PRECISION XDEBUT,XFIN,XSECT(NSMAX),XTMAIL(LMAX),DEL,DXB
+
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+ NUMSEC=0
+ JFIN=1
+
+130 NUMSEC=NUMSEC+1
+ IF(NUMSEC.LE.NBSECT) THEN
+
+
+ DXS =0.5*(DXSECT(NUMSEC)+DXSEC2(NUMSEC))
+ XDEBUT=XSECT(NUMSEC)
+ XFIN=XSECT(NUMSEC+1)
+ LMSEC=NINT((XFIN-XDEBUT)/DXS)
+ DXE = (XFIN-XDEBUT)/FLOAT(LMSEC)
+ IF(LMSEC.GT.1)THEN
+ DEL=(DXSEC2(NUMSEC)-DXSECT(NUMSEC))/FLOAT(LMSEC-1)
+ ELSEIF(LMSEC.LT.1)THEN
+ WRITE(*,*)'ERREUR DANS MAIL0,SECTEUR NUMERO ',NUMSEC
+ STOP
+ ENDIF
+ DXB=DXE-0.5*FLOAT(LMSEC-1)*DEL
+C -- CREATION DU TABLEAU DES COORDONNEES DES INTERFACES
+ JDEBUT=JFIN
+ XTMAIL(JDEBUT)=XDEBUT
+ DO 131 I=1,LMSEC
+ XTMAIL(JDEBUT+I)=XTMAIL(JDEBUT)+DXB*
+ + FLOAT(I)+0.5*DEL*FLOAT(I*(I-1))
+C XTMAIL(JDEBUT+I)=XTMAIL(JDEBUT)+(XFIN-XDEBUT)*
+C + FLOAT(I)/FLOAT(LMSEC)
+131 CONTINUE
+ JFIN=JDEBUT+LMSEC
+ GOTO 130
+ ELSE
+ LM=JFIN+1
+ ENDIF
+
+C -- REMPLISSAGE DE TMAIL
+
+ TMAIL(1)=XTMAIL(1)
+ DO 140 J=2,LM-1
+ TMAIL(J)=(XTMAIL(J)+XTMAIL(J-1))/2.
+140 CONTINUE
+ TMAIL(LM)=XTMAIL(LM-1)
+
+C DO 1 J=1,LM
+C WRITE(*,*)XTMAIL(J),TMAIL(J)
+C1 CONTINUE
+
+ END
+
+
+C--------------------------------------------------------------C
+ SUBROUTINE SAIDEV
+C--------------------------------------------------------------C
+C SON ROLE : saisie du deversement lateral C
+C--------------------------------------------------------------C
+C IMPLICIT NONE
+
+ INTEGER NSMAX,LMAX
+ PARAMETER(LMAX=3000,NSMAX=3000)
+
+ INTEGER CHOIX,ENCORE,I,J,IODEV,LM,NBDEV,N1(NSMAX),N2(NSMAX)
+ :,N3(NSMAX),IMAX,IX,N
+
+ DOUBLE PRECISION MU0(NSMAX),XDEV(NSMAX),YD0(NSMAX),EXPOS0(NSMAX)
+ DOUBLE PRECISION MU(LMAX),TMAIL(LMAX),XTMAIL(LMAX),YD(LMAX)
+ :,EXPOS(LMAX)
+ DOUBLE PRECISION QMU(NSMAX,NSMAX),TMU(NSMAX,NSMAX),X1(NSMAX)
+ :,X2(NSMAX),DMU(NSMAX),SMU(NSMAX)
+ :,QSMU(NSMAX,NSMAX),TSMU(NSMAX,NSMAX)
+ CHARACTER ETUDE*20,NOMFIC*38,CH1*1
+
+ COMMON/ETUDE/ETUDE
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+1000 WRITE(*,*)'VOULEZ VOUS CARACTERISER DU DEVERSEMENT LATERAL ?'
+C WRITE (*,*)' OUI --> 1 NON --> 0'
+C READ(*,*) IODEV
+
+C IF (IODEV.EQ.1) THEN
+ WRITE(*,*)'PAR SAISIE D''UN FICHIER dever0 : CHOIX 1'
+ WRITE(*,*)'PAR LECTURE D''UN FICHIER dever0 : CHOIX 2'
+ WRITE(*,*)'PAR SAISIE D''UN FICHIER hydev0 : CHOIX 3'
+ WRITE(*,*)'PAR LECTURE D''UN FICHIER hydev0 : CHOIX 4'
+ WRITE(*,*)'NE RIEN FAIRE : CHOIX 5'
+ READ(*,*) CHOIX
+ GOTO (1,2,3,4,5) CHOIX
+
+C 1) SAISIE PAR ECRITURE DU FICHIER DEVER0
+C ----------------------------------------
+
+C 11)SAISIE
+C----------
+C1 ENCORE = 1
+C 12)SAUVEGARDE DANS LE FICHIER DEVER0
+C-------------------------------------------
+
+C WRITE(*,*)'VOULEZ VOUS SAUVEGARDER DANS UN FICHIER ?'
+C WRITE(*,*)' OUI --> 1 NON --> 0'
+C READ (*,*) ENCORE
+C IF (ENCORE.EQ.1) THEN
+C WRITE(25,*) NBDEV
+C DO 14 I=1,NBDEV
+C WRITE(25,*)XDEV(I),YD0(I),MU0(I),EXPOS0(I)
+C14 CONTINUE
+
+ 1 I = 1
+ NOMFIC='dever0.'//ETUDE
+ OPEN(25,FILE=NOMFIC,FORM='FORMATTED',STATUS='UNKNOWN')
+C IF (ENCORE.EQ.1) THEN
+ WRITE(*,*)'LES PARAMETRES RENTRES SONT INTERPOLES'
+ WRITE(*,*)'LINEAIREMENT AUX CENTRES DE MAILLE.'
+ WRITE(*,*)'POUR UN DEVERSEMENT NUL SUR UNE ZONE'
+ WRITE(*,*)'IL FAUT RENTRER UN COEFFICIENT MU NUL'
+ WRITE(*,*)'AUX 2 EXTREMITES DE LA ZONE'
+ write(*,*)'ENTREZ DEUX FOIS MEME ABSCISSE POUR SAUT DE MU'
+ WRITE(*,*)'ABSCISSE AMONT ZONE 1 = AMONT MODELE'
+ WRITE(*,*)'ABSCISSE AVAL DERNIERE ZONE = AVAL MODELE'
+ 12 WRITE(*,*)'ENTREZ L''ABSCISSE AVAL DE LA ZONE DE DEVERSEMENT :'
+ READ (*,*)XDEV(I)
+ WRITE(*,*)'ENTREZ LA COTE DE DEVERSEMENT PUIS'
+ WRITE(*,*)'LE COEFFICIENT MU'
+ READ (*,*)YD0(I),MU0(I)
+ WRITE(*,*)'VOULEZ VOUS RENTRER UNE PUISSANCE AUTRE QUE 1.5?'
+ WRITE(*,*)' OUI --> 1 NON --> 0 OU AUTRE'
+ READ (*,'(A1)') CH1
+ IF(CH1.EQ.'1')THEN
+ WRITE(*,*)'RENTREZ L''EXPOSANT'
+ READ(*,*)EXPOS0(I)
+ ELSE
+ EXPOS0(I)=1.5
+ ENDIF
+ WRITE(25,*)XDEV(I),YD0(I),MU0(I),EXPOS0(I)
+ WRITE(*,*)'AVEZ VOUS FINI LA SAISIE DES DEVERSEMENTS ?'
+ WRITE(*,*)' OUI --> 1 NON --> 0 OU AUTRE'
+ READ (*,'(A1)') CH1
+C ENDIF
+ IF(CH1.NE.'1')THEN
+ I = I+1
+ GOTO 12
+ ENDIF
+ NBDEV = I
+ CLOSE(25)
+ GOTO 99
+
+C 2)-->SAISIE PAR LECTURE DU FICHIER DEVER0
+C ------------------------------------------------
+
+2 NOMFIC = 'dever0.'//ETUDE
+ OPEN(25,FILE=NOMFIC,FORM='FORMATTED',STATUS='OLD')
+C READ (25,*)NBDEV
+ DO 22 I=1,LMAX
+ READ(25,*,END=722)XDEV(I),YD0(I),MU0(I),EXPOS0(I)
+22 CONTINUE
+ 722 CLOSE(25)
+ NBDEV=I-1
+C 3)--> CHOIX D'UNE INTERPOLATION
+C -----------------------------
+
+99 WRITE(*,*)'VOULEZ VOUS INTERPOLER ?'
+ WRITE(*,*)' OUI --> 1 NON --> 0'
+ READ (*,*) ENCORE
+ IF (ENCORE.EQ.1) THEN
+ NOMFIC = 'devers.'//ETUDE
+ OPEN(26,FILE=NOMFIC,FORM='FORMATTED',STATUS='UNKNOWN')
+ DO 398 J=1,NBDEV-1
+ IF(XDEV(J).GT.XDEV(J+1))THEN
+ WRITE(*,*)'ABSCISSES NON CROISSANTES DANS DEVER0'
+ STOP
+ ENDIF
+398 CONTINUE
+ NOMFIC='mail.'//ETUDE
+ OPEN(20,FILE=NOMFIC,STATUS='OLD')
+ READ (20,*) LM
+ READ (20,*) (TMAIL(I),I=1,LM)
+ READ (20,*) (XTMAIL(I),I=1,LM-1)
+ CLOSE(20)
+
+ CALL INTERP(XDEV,NBDEV,YD,YD0)
+ CALL INTERP(XDEV,NBDEV,MU,MU0)
+ CALL INTERP(XDEV,NBDEV,EXPOS,EXPOS0)
+ WRITE(26,*)LM
+ DO 24 I=1,LM
+ WRITE(26,*)I,YD(I),MU(I),EXPOS(I)
+24 CONTINUE
+ CLOSE (26)
+ ENDIF
+C ENDIF
+ GO TO 1000
+C RETURN
+C SAISIE DU FICHIER HYDEV0
+ 3 OPEN(26,FILE='hydev0.'//ETUDE,STATUS='UNKNOWN')
+ I=0
+ 759 I=I+1
+ WRITE(*,*)'DONNEZ POUR LE SECTEUR NUMERO',I
+ WRITE(*,*)'LES ABSCISSES AMONT ET AVAL'
+ WRITE(*,*)'ET LE NOMBRE DE POINTS DE L''HYDROGRAMME'
+ READ(*,*)X1(I),X2(I),N3(I)
+ write(*,*)'donnez le diametre et etendue sediments'
+ read(*,*)dmu(i),smu(i)
+cliq dmu(i)=0.1
+cliq smu(i)=1.
+ WRITE(26,*)X1(I),X2(I),N3(I),dmu(i),smu(i)
+ WRITE(*,*)'DONNEZ LES COUPLES (Q,T)DANS L''ORDRE CHRONOLOGIQUE'
+ DO J=1,N3(I)
+ WRITE(*,*)J,'eme couple'
+ READ(*,*)QMU(I,J),TMU(I,J)
+ ENDDO
+ WRITE(*,*)'Donnez les couples (Qs,t) dans l''ordre chronologique'
+ DO J=1,N3(I)
+ WRITE(*,*)J,'eme couple'
+ READ(*,*)QSMU(I,J),TSMU(I,J)
+cliq qsmu(i,j)=0.
+cliq tsmu(i,j)=tmu(i,j)
+ ENDDO
+ DO J=1,N3(I)
+ WRITE(26,*)QMU(I,J),TMU(I,J),QSMU(I,J),TSMU(I,J)
+ ENDDO
+ WRITE(*,*)'VOULEZ-VOUS UN SECTEUR DE PLUS ?'
+ WRITE(*,*)'OUI --->1 NON --->0'
+ READ(*,*)IODEV
+ IF(IODEV.EQ.1)GO TO 759
+ CLOSE(26)
+ IMAX=I
+ GO TO 798
+C LECTURE DANS HYDEV0
+ 4 OPEN(26,FILE='hydev0.'//ETUDE,STATUS='OLD')
+ DO 3758 I=1,NSMAX
+ READ(26,*,END=799)X1(I),X2(I),N3(I),dmu(i),smu(i)
+ DO 4758 J=1,N3(I)
+ READ(26,*)QMU(I,J),TMU(I,J),QSMU(I,J),TSMU(I,J)
+ 4758 CONTINUE
+ 3758 CONTINUE
+ 799 CLOSE (26)
+ IMAX=I-1
+C INTERPOLATION DANS HYDEVS
+
+ 798 WRITE(*,*)'VOULEZ-VOUS INTERPOLER DANS UN FICHIER hydevs'
+ WRITE(*,*)'OUI ---> 1 NON ---> 0'
+ READ(*,*)ENCORE
+ IF(ENCORE.EQ.1)THEN
+ OPEN(26,FILE='hydevs.'//ETUDE,STATUS='UNKNOWN')
+ NOMFIC='mail.'//ETUDE
+ OPEN(20,FILE=NOMFIC,STATUS='OLD')
+ READ (20,*) LM
+ READ (20,*) (TMAIL(I),I=1,LM)
+ READ (20,*) (XTMAIL(I),I=1,LM-1)
+ CLOSE(20)
+ DO 5558 I=1,IMAX
+C CALCUL DES NUMEROS DE SECTION A PARTIR DES ABSCISSES
+ DO 5559 IX=3,LM-1
+ IF(XTMAIL(IX).GE.X1(I))THEN
+ N1(I)=IX
+ GO TO 5560
+ ENDIF
+ N1(I)=LM-1
+ 5559 CONTINUE
+ 5560 DO 5561 IX=3,LM-1
+ IF(XTMAIL(IX).GE.X2(I))THEN
+ N2(I)=IX
+ GO TO 5562
+ ENDIF
+ N2(I)=LM-1
+ 5561 CONTINUE
+ 5562 IF(N1(I).GT.N2(I))THEN
+ WRITE(*,*)'SECTEUR NUMERO',I, ' ERREUR'
+ WRITE(*,*)'ABSCISSE AMONT =',X1(I),' > ABSCISSE AVAL'
+ WRITE(*,*)'ABSCISSE AVAL =',X2(I)
+ WRITE(*,*)'ON INVERSE AMONT ET AVAL'
+ N=N1(I)
+ N1(I)=N2(I)
+ N2(I)=N
+c ELSEIF(N1(I).EQ.N2(I))THEN
+c WRITE(*,*)'SECTEUR NUMERO',I, ' ERREUR'
+c WRITE(*,*)'PAS DE CENTRE DE MAILLE DANS LE SECTEUR'
+c PAUSE
+c STOP
+ ENDIF
+ WRITE(26,*)N1(I),N2(I),N3(I),dmu(i),smu(i)
+ DO 5758 J=1,N3(I)
+C WRITE(26,*)QMU(I,J),TMU(I,J),0.,TMU(I,J)
+ WRITE(26,*)QMU(I,J),TMU(I,J),QSMU(I,J),TSMU(I,J)
+ 5758 CONTINUE
+ 5558 CONTINUE
+ CLOSE (26)
+ ENDIF
+C RETOUR AU MENU
+ GO TO 1000
+ 5 RETURN
+ END
+
+C***********************************************************************
+ SUBROUTINE SAIOUV
+C SAISIE DES CARACTERISTIQUES DES OUVRAGES
+C***********************************************************************
+C NOMBRE MAXI DE MAILLES ETC
+
+ INTEGER LMAX
+ PARAMETER(LMAX=3000)
+ INTEGER NOUMAX,NOEMAX
+ PARAMETER (NOUMAX=50,NOEMAX=10)
+ DOUBLE PRECISION XIA1,XIA2,ZORI
+ INTEGER IOUV,IDON,NREF,I,J,IREP
+ CHARACTER ETUDE*20,FICH*40,RIEN*9
+ LOGICAL TROUVE
+ INTEGER IA1(NOUMAX),IA2(NOUMAX)
+ : ,NOUV(NOUMAX),NBOUV
+ :,NBCOUP,K,NB1(noumax),nb2(noumax),nbi
+ DOUBLE PRECISION LONG,ZDEV
+ : ,COEF
+ : ,ZOUV,COEFIN,TRUP
+ : ,ZFERM,ZCOUP(NOEMAX*NOUMAX),QCOUP(NOEMAX*NOUMAX)
+ CHARACTER*1 TYPOUV,LIGNE*80
+
+ COMMON/ETUDE/ETUDE
+
+ TROUVE=.FALSE.
+
+ IDON=21
+ WRITE(*,*)'VOULEZ-VOUS SAISIR LES OUVRAGES'
+ WRITE(*,*)'EN ALLONGEANT UN FICHIER ouvrag.etude EXISTANT ---> 1'
+ WRITE(*,*)'OU CREER UN FICHIER ouvrag.etude NOUVEAU --> 0'
+ READ(*,*)IREP
+C WRITE(*,*)'DONNEZ LA DISTANCE DE PRECISION SUR OX ET OY'
+C WRITE(*,*)'QUI PERMET DE SAVOIR SI LES COORDONNEES SONT EXACTES'
+C READ(*,*)XDIST
+ IOUV=0
+
+ IF(IREP.EQ.1)THEN
+C LECTURE DES OUVRAGES
+ FICH='ouvrag.'//ETUDE
+
+ OPEN(IDON,FILE=FICH,STATUS='OLD',ERR=1)
+
+ DO 3 IOUV=1,NOUMAX
+C NREF :REFERENCE DES INTERFACES
+C NREF=-1:LE DEBIT DE L'OUVRAGE SE RAJOUTE AU DEBIT SAINT-VENANT
+C NREF=-2: LE SEUL DEBIT EST LE DEBIT DE L'OUVRAGE
+ READ(IDON,'(F11.3,1X,I4,1X,I4)',ERR=5,END=8)XIA1,NREF
+ :,NB1(iouv)
+ WRITE(*,*)'LECTURE DE L''OUVRAGE',IOUV
+ CALL TRXYIA(XIA1,IA1(IOUV),TROUVE)
+ IF(TROUVE)THEN
+ TROUVE=.FALSE.
+ ELSE
+ WRITE(*,*)'LES COORDONNÉES DE L''INTERFACE AMONT',
+ :' SONT FAUSSES'
+ STOP
+ ENDIF
+C NOUV NOMBRE D'OUVRAGES ELEMENTAIRES
+ READ (IDON,'(F11.3,1X,I4,1X,I4)',ERR=5)XIA2,NOUV(IOUV)
+ :,NB2(IOUV)
+ IF(NOUV(IOUV).GT.NOEMAX)THEN
+ WRITE(*,*)'TROP D''OUVRAGES ELEMENTAIRES : ',
+ :'NOMBRE LIMITE A',NOEMAX
+ ENDIF
+ CALL TRXYIA(XIA2,IA2(IOUV),TROUVE)
+ IF(TROUVE)THEN
+ TROUVE=.FALSE.
+ ELSE
+ WRITE(*,*)'LES COORDONNÉES DE L''INTERFACE AVAL',
+ :' NE CORRESPONDENT PAS A UNE INTERFACE'
+ STOP
+ ENDIF
+C DO 7 I=IA1(IOUV),IA2(IOUV),1
+C NREFA(I)=NREF
+C 7 CONTINUE
+C LECTURE DES CARACTERISTIQUES DES OUVRAGES
+ DO 4 J=1,NOUV(IOUV)
+C FORMAT FLUVIA
+ READ(IDON,'(A1,A9,F10.0,F10.0,F10.0,F10.0)',ERR=5)
+ :TYPOUV,RIEN,LONG,ZDEV,ZORI
+ :,COEF
+ IF(TYPOUV.EQ.'D'.OR.TYPOUV.EQ.'d')THEN
+C TYPOUV='D'
+C O correspond au meme ouvrage que D
+C mais avec une longueur et en circulaire
+ ELSEIF(TYPOUV.EQ.'O'.OR.TYPOUV.EQ.'o')THEN
+c TYPOUV(IOUV,J)='O'
+ ELSEIF(TYPOUV.EQ.'Y'.OR.TYPOUV.EQ.'y')THEN
+C TYPOUV='Y'
+ READ(IDON,*)COEFIN,ZOUV,ZFERM
+ ELSEIF(TYPOUV.EQ.'Z'.OR.TYPOUV.EQ.'z')THEN
+ READ(RIEN,'(I9)')NBCOUP
+ IF(NBCOUP.LT.1.OR.NBCOUP.GT.NOUMAX*NOEMAX)
+ :THEN
+ WRITE(*,*)'OUVRAGE DE TYPE Z'
+ WRITE(*,*)'ERREUR DANS LE NOMBRE DE COUPLES'
+ WRITE(*,*)'LE NOMBRE DE COUPLES EST: ',NBCOUP
+ WRITE(*,*)'ALORS QUE LE NOMBRE MAXIMAL PERMIS EST :',
+ :NOUMAX*NOEMAX
+ STOP
+ ENDIF
+ DO 9 K=1,NBCOUP
+ READ(IDON,*,ERR=5)ZCOUP(K),QCOUP(K)
+ 9 CONTINUE
+ ELSEIF(TYPOUV.EQ.'B'.OR.TYPOUV.EQ.'b')THEN
+ READ(IDON,'(4F10.2)',ERR=5) LONG,ZDEV,ZORI,COEF
+ READ(IDON,'(2F10.2)',ERR=5)LONG,ZDEV
+ READ(IDON,'(2F6.2,I2,I4)',ERR=5) LONG,ZDEV,K,NBCOUP
+c ELSEIF(TYPOUV.EQ.'T'.OR.TYPOUV.EQ.'t')THEN
+c READ(IDON,'(4F10.2)',ERR=5) LONG,ZDEV,ZORI,COEF
+c READ(IDON,'(2F10.2)',ERR=5)LONG,ZDEV
+c READ(IDON,'(2F6.2,I2,I4)',ERR=5) LONG,ZDEV,K,NBCOUP
+ ELSE
+ WRITE(*,*)'ERREUR DANS LE TYPE D''OUVRAGE'
+ ENDIF
+C HAUT(IOUV,J)=ZORI-ZDEV(IOUV,J)
+ 4 CONTINUE
+ 3 CONTINUE
+ WRITE(*,*)'ATTENTION'
+ WRITE(*,*)'LE NOMBRE MAXIMAL D''OUVRAGES ',NOUMAX,' A ETE LU'
+ 8 NBOUV=IOUV-1
+ GO TO 6
+ 1 WRITE(*,*)'LE FICHIER ouvrag.',ETUDE,' N''EXISTE PAS'
+ NBOUV=0
+ GO TO 6
+ 5 WRITE (*,*)'ERREUR DANS LA LECTURE DU FICHIER ouvrag.',ETUDE
+ STOP
+ 6 CONTINUE
+ WRITE(*,*)'VOULEZ-VOUS RENTREZ UN OUVRAGE SUPPLEMENTAIRE'
+ WRITE(*,*)'OUI ---> 1 NON --->0'
+ READ(*,*)IREP
+ IF(IREP.EQ.1)THEN
+C ON VA A LA SAISIE
+ GO TO 120
+ ENDIF
+C ON VA EN FIN DE PROGRAMME
+ GO TO 110
+C SAISIE DES OUVRAGES
+ ELSEIF(IREP.EQ.0)THEN
+ FICH='ouvrag.'//ETUDE
+ OPEN(IDON,FILE=FICH,STATUS='UNKNOWN')
+ ENDIF
+ 120 IOUV=IOUV+1
+ IF(IOUV.GE.NOUMAX)THEN
+ GO TO 110
+ ELSE
+ WRITE(*,*)'OUVRAGE ',IOUV
+ 130 WRITE(*,*)'DONNEZ LES COORDONNEES DE L''ABSCISSE AMONT'
+ WRITE(*,*)'XIA1 ?'
+ READ(*,*)XIA1
+ CALL TRXYIA(XIA1,IA1(IOUV),TROUVE)
+ IF(TROUVE)THEN
+ TROUVE=.FALSE.
+ ELSE
+ WRITE(*,*)'L''ABSCISSE AMONT EST FAUSSE'
+ GO TO 130
+ ENDIF
+ 1130 WRITE(*,*)'DONNEZ LES COORDONNEES DE L''ABSCISSE AVAL'
+ WRITE(*,*)'XIA2 ?'
+ READ(*,*)XIA2
+ CALL TRXYIA(XIA2,IA2(IOUV),TROUVE)
+ IF(TROUVE)THEN
+ TROUVE=.FALSE.
+ IF(XIA1.GT.XIA2)THEN
+ WRITE(*,*)'L''ABSCISSE AVAL EST INFERIEURE A L''AMONT'
+ GO TO 1130
+ ENDIF
+ ELSEIF(ABS(XIA2-999999.999).LT.0.0001)THEN
+ WRITE(*,*)'L''OUVRAGE DONNE A L''EXTERIEUR'
+ ELSE
+ WRITE(*,*)'L''ABSCISSE AVAL EST FAUSSE'
+ GO TO 1130
+ ENDIF
+ 102 WRITE(*,*)'DONNEZ LA REFERENCE DE L''OUVRAGE : -1 OU -2'
+ WRITE(*,*)'-1 SIGNIFIE LE DEBIT S''AJOUTE A SAINT-VENANT'
+ WRITE(*,*)'-2 SIGNIFIE : SEUL LE DEBIT DE L''OUVRAGE EST PRIS'
+ READ(*,*)NREF
+ IF(NREF.NE.-1.AND.NREF.NE.-2) GO TO 102
+ write(*,*)'donnez le numero de bief abscisse amont',
+ :' (touche entree si sans objet)'
+ read(*,*)ligne
+ If (ligne(1:1).eq.' ')then
+ WRITE(IDON,'(F11.3,1X,I4,A5)')XIA1,NREF,' '
+ else
+ read(ligne,*)nbi
+ WRITE(IDON,'(F11.3,1X,I4,1X,I4)')XIA1,NREF,nbi
+ endif
+
+ WRITE(*,*)'DONNEZ LE NOMBRE D''OUVRAGES ELEMENTAIRES'
+ READ(*,*)NOUV(IOUV)
+ if(nouv(iouv).gt.noemax)then
+ write(*,*)'nombre ouvrages elementaires limite a ',noemax
+ nouv(iouv)=10
+ endif
+ write(*,*)'donnez le numero de bief abscisse aval',
+ :' (touche entree si sans objet)'
+ read(*,*)ligne
+ If (ligne(1:1).eq.' ')then
+ WRITE(IDON,'(F11.3,1X,I4,A5)')XIA2,NOUV(IOUV),' '
+ else
+ read(ligne,*)nbi
+ WRITE(IDON,'(F11.3,1X,I4,1X,I4)')XIA2,NOUV(IOUV),nbi
+ endif
+C WRITE(IDON,*)XIA2,NOUV(IOUV)
+ DO 108 I=1,NOUV(IOUV)
+ WRITE(*,*)'DONNEZ LE TYPE D''OUVRAGE'
+ WRITE(*,*) 'B : SIGNIFIE DIGUE ERODABLE BRECHE '
+ WRITE(*,*) 'D : SIGNIFIE DEVERSOIR ORIFICE'
+ WRITE(*,*) 'O : SIGNIFIE ORIFICE CRCULAIRE'
+ WRITE(*,*) 'Y : SIGNIFIE DEVERSOIR ORIFICE MOBILE'
+ WRITE(*,*) 'Z : SIGNIFIE LOI DEBIT COTE DONNEE'
+c WRITE(*,*) 'T : SIGNIFIE DIGUE ERODABLE BRECHE TRAPEZE'
+
+ 109 READ(*,'(A1)')TYPOUV
+ IF (TYPOUV.NE.'D'.AND.TYPOUV.NE.'d'.AND.
+ : TYPOUV.NE.'Z'.AND.TYPOUV.NE.'z'.AND.
+ : TYPOUV.NE.'Q'.AND.TYPOUV.NE.'q'.AND.
+ : TYPOUV.NE.'B'.AND.TYPOUV.NE.'b'.AND.
+ : TYPOUV.NE.'O'.AND.TYPOUV.NE.'o'.AND.
+c : TYPOUV.NE.'T'.AND.TYPOUV.NE.'T'.AND.
+ : TYPOUV.NE.'Y'.AND.TYPOUV.NE.'y')THEN
+ WRITE(*,*)'OUVRAGE INCONNU OU NON DISPONIBLE'
+ GO TO 109
+ ENDIF
+ IF(TYPOUV.EQ.'Y'.OR.TYPOUV.EQ.'y')THEN
+ TYPOUV='Y'
+ WRITE(*,*)'DONNEZ LA LONGUEUR DE SEUIL'
+ READ(*,*)LONG
+ WRITE(*,*)'DONNEZ LA COTE DE SEUIL'
+ READ(*,*)ZDEV
+ WRITE(*,*)'DONNEZ LA COTE DE MISE EN CHARGE'
+ READ(*,*)ZORI
+ WRITE(*,*)'DONNEZ LE COEFFICIENT DE DEBIT'
+ READ(*,*)COEF
+ WRITE(IDON,'(A1,9X,4F10.3)')TYPOUV,LONG,ZDEV,ZORI,COEF
+ WRITE(*,*)'DONNEZ LE COEFFICIENT INITIAL'
+ READ(*,*)COEFIN
+ WRITE(*,*)'DONNEZ LA COTE D''OUVERTURE'
+ READ(*,*)ZOUV
+ WRITE(*,*)'DONNEZ LA COTE DE FERMETURE'
+ READ(*,*)ZFERM
+ WRITE(IDON,*)COEFIN,ZOUV,ZFERM
+ ELSEIF(TYPOUV.EQ.'D'.OR.TYPOUV.EQ.'d')THEN
+ TYPOUV='D'
+ WRITE(*,*)'DONNEZ LA LONGUEUR DE SEUIL'
+ READ(*,*)LONG
+ WRITE(*,*)'DONNEZ LA COTE DE SEUIL'
+ READ(*,*)ZDEV
+ WRITE(*,*)'DONNEZ LA COTE DE MISE EN CHARGE'
+ READ(*,*)ZORI
+ WRITE(*,*)'DONNEZ LE COEFFICIENT DE DEBIT'
+ READ(*,*)COEF
+ WRITE(IDON,'(A1,9X,4F10.3)')TYPOUV,LONG,ZDEV,ZORI,COEF
+C O correspond au meme ouvrage que D mais
+C avec une longueur et en circulaire
+ ELSEIF(TYPOUV.EQ.'O'.OR.TYPOUV.EQ.'o')THEN
+ TYPOUV='O'
+ WRITE(*,*)'DONNEZ LA LONGUEUR DE CONDUITE'
+ READ(*,*)LONG
+ WRITE(*,*)'DONNEZ LA COTE INFERIEURE'
+ READ(*,*)ZDEV
+ WRITE(*,*)'DONNEZ LE DIAMETRE DE LA CONDUITE'
+ READ(*,*)ZORI
+ WRITE(*,*)'DONNEZ LE COEFFICIENT DE DEBIT'
+ READ(*,*)COEF
+ WRITE(IDON,'(A1,9X,4F10.3)')TYPOUV,LONG,ZDEV,ZORI,COEF
+ ELSEIF(TYPOUV.EQ.'Z'.OR.TYPOUV.EQ.'z')THEN
+ TYPOUV='Z'
+ 7000 WRITE(*,*)'DONNEZ LE NOMBRE DE COUPLES'
+ READ(*,*)NBCOUP
+ IF(NBCOUP.GT.NOEMAX*NOUMAX)THEN
+ WRITE(*,*)'TROP DE COUPLES'
+ WRITE(*,*)'DONNEZ UN NOMBRE PLUS FAIBLE'
+ GO TO 7000
+ ENDIF
+ WRITE(IDON,'(A1,I9,4F10.3)')TYPOUV,NBCOUP,0.,0.,0.,0.
+ WRITE(*,*)'ATTENTION : LES COTES DOIVENT AUGMENTER'
+ DO K=1,NBCOUP
+ WRITE(*,*)'DONNEZ LE COUPLE(Z,Q)',K
+ READ(*,*)ZCOUP(K),QCOUP(K)
+ WRITE(IDON,'(2F10.3)')ZCOUP(K),QCOUP(K)
+ ENDDO
+ ELSEIF(TYPOUV.EQ.'Q'.OR.TYPOUV.EQ.'q')THEN
+ TYPOUV='Q'
+ 7001 WRITE(*,*)'DONNEZ LE NOMBRE DE COUPLES'
+ READ(*,*)NBCOUP
+ IF(NBCOUP.GT.NOEMAX*NOUMAX)THEN
+ WRITE(*,*)'TROP DE COUPLES'
+ WRITE(*,*)'DONNEZ UN NOMBRE PLUS FAIBLE'
+ GO TO 7001
+ ENDIF
+ WRITE(IDON,'(A1,I9,4F10.3)')TYPOUV,NBCOUP,0.,0.,0.,0.
+ WRITE(*,*)'ATTENTION : LES COTES DOIVENT AUGMENTER'
+ DO K=1,NBCOUP
+ WRITE(*,*)'DONNEZ LE COUPLE(T,Q)',K
+ READ(*,*)ZCOUP(K),QCOUP(K)
+ WRITE(IDON,'(2F10.3)')ZCOUP(K),QCOUP(K)
+ ENDDO
+ ELSEIF(TYPOUV.EQ.'B'.OR.TYPOUV.EQ.'b')THEN
+ TYPOUV='B'
+ WRITE(*,*)'DONNEZ LA COTE DE CRETE'
+ READ(*,*)LONG
+ WRITE(*,*)'DONNEZ LE temps de debut de rupture'
+ READ(*,*)TRUP
+ WRITE(*,*)'DONNEZ LA COTE DE PIED'
+ READ(*,*)ZDEV
+ WRITE(*,*)'DONNEZ LA LARGEUR EN CRETE'
+ READ(*,*)ZORI
+ WRITE(*,*)'DONNEZ LA LARGEUR EN PIED'
+ READ(*,*)COEF
+ WRITE(IDON,'(A1,1x,F8.0,4F10.3)')
+ :TYPOUV,TRUP,LONG,ZDEV,ZORI,COEF
+ WRITE(*,*)'DONNEZ LE D50 (en mm)'
+ READ(*,*)LONG
+ WRITE(*,*)'DONNEZ LE STRICKLER DE BRECHE'
+ READ(*,*)ZDEV
+ WRITE(*,*)'DONNEZ LA MASSE VOLUMIQUE (en KG/M3)'
+ READ(*,*)ZORI
+ WRITE(*,*)'DONNEZ LA POROSITE (<1)'
+ READ(*,*)COEF
+ WRITE(IDON,'(4F10.2)')LONG,ZDEV,ZORI,COEF
+ WRITE(*,*)'DONNEZ LA COTE DE FOND DE BRECHE INITIALE'
+ READ(*,*)LONG
+ WRITE(*,*)'DONNEZ LA LARGEUR DE BRECHE INITIALE(mm)'
+ WRITE(*,*)'RAYON OU LARGEUR'
+ READ(*,*)ZDEV
+ WRITE(*,*)'DONNEZ LA LARGEUR DE BRECHE MAXIMALE(m)'
+ READ(*,*)ZORI
+ WRITE(IDON,'(3F10.2)')LONG,ZDEV,ZORI
+ WRITE(*,*)'DONNEZ LE PAS DE TEMPS DE STOCKAGE'
+ WRITE(*,*)'DES RESULTATS SUR LA BRECHE'
+ READ(*,*)LONG
+ WRITE(*,*)'DONNEZ LE COEFFICIENT DE PERTE DE CHARGE'
+ WRITE(*,*)'AMONT DE LA BRECHE (<0.5)'
+ READ(*,*)ZDEV
+ WRITE(*,*)'DONNEZ LE TYPE DE RUPTURE'
+ WRITE(*,*)'RENARD ---> 0'
+ WRITE(*,*)'SURVERSE--->-1'
+ READ(*,*)K
+ WRITE(*,*)'DONNEZ LE NOMBRE MAXIMAL DE PAS DE TEMPS'
+ WRITE(*,*)'DES RESULTATS SUR LA BRECHE (<9000)'
+ READ(*,*)NBCOUP
+ WRITE(IDON,'(2F6.2,I2,I4)') LONG,ZDEV,K,NBCOUP
+c ELSEIF(TYPOUV.EQ.'T'.OR.TYPOUV.EQ.'t')THEN
+c TYPOUV='T'
+c WRITE(*,*)'DONNEZ LA COTE DE CRETE'
+c READ(*,*)LONG
+c WRITE(*,*)'DONNEZ LA COTE DE PIED'
+c READ(*,*)ZDEV
+c WRITE(*,*)'DONNEZ LA LARGEUR EN CRETE'
+c READ(*,*)ZORI
+c WRITE(*,*)'DONNEZ LA LARGEUR EN PIED'
+c READ(*,*)COEF
+c WRITE(IDON,'(A1,9X,4F10.3)')TYPOUV,LONG,ZDEV,ZORI,COEF
+c WRITE(*,*)'DONNEZ LE D50 (en mm)'
+c READ(*,*)LONG
+c WRITE(*,*)'DONNEZ LE STRICKLER DE BRECHE'
+c READ(*,*)ZDEV
+c WRITE(*,*)'DONNEZ LA MASSE VOLUMIQUE (en KG/M3)'
+c READ(*,*)ZORI
+c WRITE(*,*)'DONNEZ LA POROSITE (<1)'
+c READ(*,*)COEF
+c WRITE(IDON,'(4F10.2)')LONG,ZDEV,ZORI,COEF
+c WRITE(*,*)'DONNEZ LA COTE DE FOND DE BRECHE INITIALE'
+c READ(*,*)LONG
+c WRITE(*,*)'DONNEZ LA LARGEUR DE BRECHE INITIALE(mm)'
+c WRITE(*,*)'RAYON OU LARGEUR'
+c READ(*,*)ZDEV
+c WRITE(IDON,'(2F10.2)')LONG,ZDEV
+c WRITE(*,*)'DONNEZ LE PAS DE TEMPS DE STOCKAGE'
+c WRITE(*,*)'DES RESULTATS SUR LA BRECHE'
+c READ(*,*)LONG
+c WRITE(*,*)'DONNEZ LE COEFFICIENT DE PERTE DE CHARGE'
+c WRITE(*,*)'AMONT DE LA BRECHE (<0.5)'
+c READ(*,*)ZDEV
+c WRITE(*,*)'DONNEZ LE TYPE DE RUPTURE'
+c WRITE(*,*)'RENARD ---> 0'
+c WRITE(*,*)'SURVERSE--->-1'
+c READ(*,*)K
+c WRITE(*,*)'DONNEZ LE NOMBRE MAXIMAL DE PAS DE TEMPS'
+c WRITE(*,*)'DES RESULTATS SUR LA BRECHE (<1999)'
+c READ(*,*)NBCOUP
+c WRITE(IDON,'(2F6.2,I2,I4)') LONG,ZDEV,K,NBCOUP
+C fin du if sur le type d'ouvrage
+ ENDIF
+ 108 CONTINUE
+ WRITE(*,*)'VOULEZ-VOUS RENTREZ UN NOUVEL OUVRAGE'
+ WRITE(*,*)'OUI ---> 1 NON --->0'
+ READ(*,*)IREP
+ IF(IREP.EQ.1)THEN
+ GO TO 120
+ ENDIF
+C FIN DU IF SUR NOUMAX DEPASSE
+ ENDIF
+ 110 CLOSE(IDON)
+ RETURN
+ END
+
+
+C------------------------------------------------------------------C
+ SUBROUTINE SAISAV
+C------------------------------------------------------------------C
+C SON ROLE : ecriture des fichiers contenant les parametres C
+C de sauvegarde C
+C------------------------------------------------------------------C
+C maille ---> abscisses de sauvegarde des hydrogrammes C
+C limnigrammes C
+C vitigrammes C
+C tsor -----> temps de sauvegarde des profils de debit C
+C hauteur C
+C vitesse C
+C------------------------------------------------------------------C
+
+C IMPLICIT NONE
+
+ INTEGER J,NBSSAV,NTSOR
+ DOUBLE PRECISION MAILLE(100),TSOR(100)
+
+ CHARACTER *20 ETUDE,NOMFIC*38
+
+ COMMON/ETUDE/ETUDE
+
+
+ WRITE(*,*)'ENTREZ LE NOMBRE D''HYDROGRAMMES DESIRES(MAX=100)'
+ READ(*,*) NBSSAV
+ IF ((NBSSAV.GT.0).AND.(NBSSAV.LE.100)) THEN
+ NOMFIC='abshyd.'//ETUDE
+ OPEN(75,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(75,*) NBSSAV
+ DO 10 J=1,NBSSAV
+ WRITE(*,*)J,'IEME ABSCISSE :'
+ READ(*,*) MAILLE(J)
+ WRITE(75,'(I4,1X,F11.3,1X,I4)')J,MAILLE(J),1
+10 CONTINUE
+ CLOSE(75)
+ ENDIF
+
+ WRITE(*,*)'ENTREZ LE NOMBRE DE PROFILS DESIRES(MAX=100)'
+ READ(*,*) NTSOR
+ IF ((NTSOR.GT.0).AND.(NTSOR.LE.100)) THEN
+ NOMFIC='tnprof.'//ETUDE
+ OPEN(76,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(76,*)NTSOR
+ DO 20 J=1,NTSOR
+ WRITE(*,*)J,'IEME TEMPS DE SAUVEGARDE :'
+ READ(*,*)TSOR(J)
+ WRITE(76,*)J,TSOR(J)
+20 CONTINUE
+ CLOSE(76)
+ ENDIF
+
+ RETURN
+ END
+
+
+C**********************************************************************
+ SUBROUTINE TRXYIA (X,I,TROUVE)
+C---------------------------------------------------------------------C
+C VERIFIE QUE LES COORDONNEES DES INTERFACES D'UN OUVRAGE SONT BONNES
+C ENTREE: X DE L'INTERFACE
+C SORTIE: NUMERO DE L'INTERFACE
+C LA VARIABLE LOGIQUE TROUVE INDIQUE SI L'INTERFACE EST TROUVEE
+C=====================================================================C
+ INTEGER LMAX,LM
+ PARAMETER (LMAX=3000)
+ DOUBLE PRECISION TMAIL(LMAX),XTMAIL(LMAX)
+ CHARACTER ETUDE*20,NOMFIC*38
+ INTEGER I,N
+ DOUBLE PRECISION X
+ LOGICAL TROUVE
+
+ COMMON/ETUDE/ETUDE
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+C -----------------LECTURE DU FICHIER MAIL------------------------C
+
+ N=20
+ NOMFIC='mail.'//ETUDE
+ OPEN(N,FILE=NOMFIC,STATUS='OLD')
+ READ (N,*) LM
+ READ (N,*) (TMAIL(I),I=1,LM)
+ READ (N,*) (XTMAIL(I),I=1,LM-1)
+ CLOSE(N)
+
+C recherche du xtmail correspondant
+ DO 1 I=1,LM-1
+C 0.006 CAR LES COORDONNEES SONT PRECISES AU CM
+ IF(ABS(X-XTMAIL(I)).LT.0.006)THEN
+ TROUVE=.TRUE.
+ GO TO 2
+ ENDIF
+ 1 CONTINUE
+ 2 RETURN
+ END
+C-------------------------------------------------
+ DOUBLE PRECISION FUNCTION DETL(Y,I)
+C--------------------------------------------------
+C DETERMINATION DE LA LARGEUR L(Y)
+C--------------------------------------------------
+C IMPLICIT NONE
+
+ INTEGER NSMAX,NCMAX,I,N,IJ,IE
+
+ PARAMETER(NSMAX=3000,NCMAX=300)
+
+ INTEGER NC0(NSMAX),IDPARS
+
+ CHARACTER*20 ETUDE
+
+ DOUBLE PRECISION Y,EPS,YI,LI,
+ + LISEC0(NCMAX,NSMAX),YISEC0(NCMAX,NSMAX),CTDF0(NSMAX)
+
+ COMMON/NC0/NC0
+ COMMON/REEL/YISEC0,LISEC0
+ COMMON/CTD/CTDF0
+ COMMON/ETUDE/ETUDE
+ COMMON/CONST/IDPARS,IE
+ COMMON/EPSI/EPS
+
+ DO 2 N=2,NC0(I)
+ IF (Y.LT.YISEC0(N,I)) THEN
+ IJ = N-1
+ GO TO 3
+ ENDIF
+2 CONTINUE
+ IJ = NC0(I)
+3 YI=YISEC0(IJ,I)
+ LI=LISEC0(IJ,I)
+ IF((Y-YI).GT.EPS)THEN
+ DETL = LI+(Y-YI)*(LISEC0(IJ+1,I)-LI)/(YISEC0(IJ+1,I)-YI)
+ ELSE
+ DETL=LI
+ ENDIF
+ RETURN
+ END
+C-----------------------------------------------------------------------
+ SUBROUTINE GEOMET
+C-----------------------------------------------------------------------
+C Menu relatif à la saisie de la géométrie
+C-----------------------------------------------------------------------
+C geomac-i00: donnees de topographie + marques de lignes directrices
+Csi elles existent
+C geomac-i0 : topographie aux intermailles +
+Cmarques de lignes directrices si elles existent
+C geomac-i : topographie aux intermailles
+C + marques directrices + description sédimentaire
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER CHOIX
+
+C Affichage du menu de saisie de la géométrie
+C-----------------------------------------------------------------------
+30 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' MENU DE SAISIE DE LA GEOMETRIE '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)'Choix 1: pour les topographies simples, saisie de '
+ WRITE(*,*)' la geometrie en abscisse-cote puis '
+ WRITE(*,*)' interpolation lineaire dans geomac-i0 '
+cliq WRITE(*,*)' interpolation lineaire dans geomac-i '
+ WRITE(*,*)' (le maillage doit pre-exister) '
+ WRITE(*,*)'Choix 2: pour les topographies plus complexes, '
+ WRITE(*,*)' ecriture de geomac-i0 a partir d''un '
+cliq WRITE(*,*)' ecriture de geomac-i a partir d''un '
+ WRITE(*,*)' fichier topographique au format Mocahy '
+ WRITE(*,*)'Choix 3: saisie des parametres sedimentaires pour '
+ WRITE(*,*)' complement de geomac-i0 en geomac-i '
+ WRITE(*,*)'Choix 4: saisie des largeurs actives par zones '
+ WRITE(*,*)'Choix 0: retour au menu principal '
+ WRITE(*,*)'--------------------------------------------------'
+ WRITE(*,*)' '
+ WRITE(*,*)' '
+ WRITE(*,*)'Votre choix:'
+ READ(*,*)CHOIX
+
+C Ventilation
+C-----------------------------------------------------------------------
+ IF(CHOIX.EQ.1)THEN
+ CALL SAITOPO
+ GOTO 30
+ ELSEIF(CHOIX.EQ.2)THEN
+ CALL TRMOCRUBE
+ GOTO 30
+ ELSEIF(CHOIX.EQ.3)THEN
+ CALL SAISED
+ RETURN
+ ELSEIF(CHOIX.EQ.4)THEN
+ CALL SAILACT
+ RETURN
+ ELSEIF(CHOIX.EQ.0)THEN
+ RETURN
+ ENDIF
+ WRITE(*,*) 'MAUVAIS CHOIX, RECOMMENCEZ'
+ GOTO 30
+ END
+C-----------------------------------------------------------------------
+ SUBROUTINE SAITOPO
+C-----------------------------------------------------------------------
+C Saisie abscisse-cote de la topographie (hors paramètres sédimentaires)
+C Écriture dans geomac-i00
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX,NSMAX,NCMAX
+ PARAMETER (LMAX=3000,NCMAX=300,NSMAX=3000)
+ INTEGER I,LM
+ INTEGER NUMCOU,NUMSEC,DERCOU
+ INTEGER NCOU(0:NSMAX)
+ DOUBLE PRECISION TMAIL(LMAX),XTMAIL(LMAX),XSECT(NSMAX)
+ DOUBLE PRECISION YCOU(NCMAX*NSMAX),ZCOU(NCMAX*NSMAX)
+ & ,YPERMUT,ZPERMUT
+ CHARACTER*1 CCOU(NCMAX*NSMAX),CPERMUT,INDIC
+ CHARACTER ETUDE*20,NOMFIC*40,REP*1
+ LOGICAL EXISTG,EXISTD
+ LOGICAL ENCORE,ENCOR2,ENCOR3
+
+ COMMON/ETUDE/ETUDE
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+C Vérification: lecture du fichier mail
+C-----------------------------------------------------------------------
+ NOMFIC='mail.'//ETUDE
+ OPEN(20,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ READ(20,*)LM
+ READ(20,*)(TMAIL(I),I=1,LM)
+ READ(20,*)(XTMAIL(I),I=1,LM-1)
+ CLOSE(20)
+
+C Preambule
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'SAISIE DE LA TOPOGRAPHIE EN ABSCISSE-COTE --->1'
+ WRITE(*,*)'ou lecture fichier geomac-i00.etude --->0'
+ Write(*,*)' entrez 0 ou 1'
+ Read (*,*)indic
+ IF(indic.eq.'0')then
+ CALL INTLINTOPO
+ RETURN
+ ENDIF
+ WRITE(*,*)' '
+ WRITE(*,*)'ATTENTION1: dans l''etape d''interpolation lineaire, '
+ & ,'les sections seront interpolees point a point; il faut donc '
+ & ,'au moins qu''il y ait le meme nombre de points par section.'
+ WRITE(*,*)' '
+ WRITE(*,*)'ATTENTION2: les section doivent être entrées de'
+ & ,' l''amont vers l''aval et avec des abscisses croissantes'
+ & ,' (eventuellement negatives).'
+ WRITE(*,*)' '
+ NOMFIC='geomac-i00.'//ETUDE
+ OPEN(24,FILE=NOMFIC,STATUS='UNKNOWN',FORM='FORMATTED')
+ NUMSEC=0
+ NCOU(0)=0
+ ENCORE=.TRUE.
+
+C Abscisse de la section
+C-----------------------------------------------------------------------
+ 2003 IF(ENCORE)THEN
+ NUMSEC=NUMSEC+1
+ 2004 WRITE(*,*)' '
+ WRITE(*,*)'On est dans la section numero ',NUMSEC
+ WRITE(*,*)'Donnez son abscisse XSECT(',NUMSEC,') = ?'
+ READ(*,*) XSECT(NUMSEC)
+ IF(NUMSEC.GE.2.AND.XSECT(NUMSEC).LE.XSECT(NUMSEC-1))THEN
+ WRITE(*,*)'ABSCISSE NON CROISSANTE, RECOMMENCEZ'
+ GOTO 2004
+ ENDIF
+ IF((XSECT(NUMSEC).LT.TMAIL(1)).OR.
+ & (XSECT(NUMSEC).GT.TMAIL(LM)))THEN
+ WRITE(*,*)'ATTENTION, section en dehors du domaine maille'
+ ENDIF
+
+C Entree des donnees geometriques pour la section
+C-----------------------------------------------------------------------
+ NUMCOU=NCOU(NUMSEC-1)
+ ENCOR2=.TRUE.
+ ENCOR3=.FALSE.
+ EXISTG=.FALSE.
+ EXISTD=.FALSE.
+ WRITE(*,*)'ATTENTION: les points abscisse-cote seront reordonnes '
+ & ,'selon des abscisses croissantes de la rive GAUCHE vers la '
+ & ,'rive DROITE (eventuellement negatives).'
+ WRITE(*,*)' '
+ 222 IF(ENCOR2)THEN
+ NUMCOU=NUMCOU+1
+ 223 WRITE(*,*)'Abscisse du ',NUMCOU-NCOU(NUMSEC-1),'ieme point'
+ & ,' dans la section = ?'
+ READ(*,*) YCOU(NUMCOU)
+ WRITE(*,*)'Cote du ',NUMCOU-NCOU(NUMSEC-1),'ieme point'
+ & ,' dans la section = ?'
+ READ(*,*) ZCOU(NUMCOU)
+ CCOU(NUMCOU)=' '
+ IF(.NOT.EXISTG)THEN
+ WRITE(*,*)'Ce point est-il la limite mineur/Moyen en rive'
+ & ,' gauche ? (O/N)'
+ READ(*,*) REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ CCOU(NUMCOU)='G'
+ EXISTG=.TRUE.
+ ENDIF
+ ELSEIF(.NOT.EXISTD)THEN
+ WRITE(*,*)'Ce point est-il la limite mineur/Moyen en rive'
+ & ,' droite ? (O/N)'
+ READ(*,*) REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ CCOU(NUMCOU)='D'
+ EXISTD=.TRUE.
+ ENDIF
+ ENDIF
+ IF(NUMCOU-NCOU(NUMSEC-1).EQ.NCMAX)THEN
+ WRITE(*,*)'C''etait le dernier point possible.'
+ WRITE(*,*)'Si la section n''est pas complete, '
+ & ,'recommencez en limitant le nombre de points a ',NCMAX
+ ENCOR2=.FALSE.
+ GOTO 222
+ ENDIF
+ IF(NUMCOU-NCOU(NUMSEC-1).GT.3)THEN
+ WRITE(*,*)'Y a-t-il d''autres points ? (O/N)'
+ READ(*,*) REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ ENCOR2=.TRUE.
+ ELSE
+ ENCOR2=.FALSE.
+ ENDIF
+ ENDIF
+ GOTO 222
+ ENDIF
+ NCOU(NUMSEC)=NUMCOU
+
+C Modification éventuelle de la section
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ 224 WRITE(*,*)'Voulez-vous modifier des points ? (O/N)'
+ READ(*,*) REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ ENCOR3=.TRUE.
+ ELSE
+ ENCOR3=.FALSE.
+ ENDIF
+ IF(ENCOR3)THEN
+ WRITE(*,*)'Numero du point a modifier ?'
+ READ(*,*)NUMCOU
+ NUMCOU=NUMCOU+NCOU(NUMSEC-1)
+ IF(NUMCOU.GE.NCOU(NUMSEC-1)+1.AND.NUMCOU.LE.NCOU(NUMSEC))THEN
+ WRITE(*,*)'Abscisse du ',NUMCOU-NCOU(NUMSEC-1),'ieme point'
+ & ,' dans la section = ?'
+ READ(*,*) YCOU(NUMCOU)
+ WRITE(*,*)'Cote du ',NUMCOU-NCOU(NUMSEC-1),'ieme point'
+ & ,' dans la section = ?'
+ READ(*,*) ZCOU(NUMCOU)
+ ELSE
+ WRITE(*,*)'MAUVAIS CHOIX, RECOMMENCEZ'
+ ENDIF
+ GOTO 224
+ ENDIF
+
+C Réarrangement des couples par abscisse croissante dans la section
+C-----------------------------------------------------------------------
+ DO DERCOU=NCOU(NUMSEC)-1,NCOU(NUMSEC-1)+1,-1
+ DO NUMCOU=NCOU(NUMSEC-1)+1,DERCOU
+ IF(YCOU(NUMCOU).GT.YCOU(NUMCOU+1))THEN
+ YPERMUT=YCOU(NUMCOU+1)
+ ZPERMUT=ZCOU(NUMCOU+1)
+ CPERMUT=CCOU(NUMCOU+1)
+ YCOU(NUMCOU+1)=YCOU(NUMCOU)
+ ZCOU(NUMCOU+1)=ZCOU(NUMCOU)
+ CCOU(NUMCOU+1)=CCOU(NUMCOU)
+ YCOU(NUMCOU)=YPERMUT
+ ZCOU(NUMCOU)=ZPERMUT
+ CCOU(NUMCOU)=CPERMUT
+ ENDIF
+ ENDDO
+ ENDDO
+
+C Ajout des limites mineur/Moyen si absentes
+C-----------------------------------------------------------------------
+ EXISTG=.FALSE.
+ EXISTD=.FALSE.
+ DO NUMCOU=NCOU(NUMSEC-1)+1,NCOU(NUMSEC)
+ IF((.NOT.EXISTG).AND.(CCOU(NUMCOU).EQ.'G'))THEN
+ EXISTG=.TRUE.
+ ELSEIF((.NOT.EXISTD).AND.(CCOU(NUMCOU).EQ.'D'))THEN
+ EXISTD=.TRUE.
+ ELSE
+ CCOU(NUMCOU)=' '
+ ENDIF
+ ENDDO
+ IF(.NOT.EXISTG) CCOU(NCOU(NUMSEC-1)+1)='G'
+ IF(.NOT.EXISTD) CCOU(NCOU(NUMSEC))='D'
+
+C Sauvegarde de la section entrée
+C-----------------------------------------------------------------------
+ WRITE(24,9992) NUMSEC,XSECT(NUMSEC),NCOU(NUMSEC)-NCOU(NUMSEC-1)
+ DO NUMCOU=NCOU(NUMSEC-1)+1,NCOU(NUMSEC)
+ WRITE(24,9991) CCOU(NUMCOU),YCOU(NUMCOU),ZCOU(NUMCOU)
+ ENDDO
+
+C Passage à la section suivante
+C-----------------------------------------------------------------------
+ IF(NUMSEC.GT.1)THEN
+ WRITE(*,*)' '
+ WRITE(*,*)'Y a-t-il d''autres sections ? (O/N)'
+ READ(*,*) REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ ENCORE=.TRUE.
+ ELSE
+ ENCORE=.FALSE.
+ ENDIF
+ ENDIF
+ GOTO 2003
+ ENDIF
+
+C Fin et sortie
+C-----------------------------------------------------------------------
+ CLOSE(24)
+ CALL INTLINTOPO
+ RETURN
+
+C Traitement des erreurs
+C-----------------------------------------------------------------------
+ 1000 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'ERREUR sur le fichier ',NOMFIC
+ WRITE(*,*)'Le fichier n''existe pas ou est illisible'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+ RETURN
+
+ 9992 FORMAT(I4,1X,F11.3,1X,I2)
+ 9991 FORMAT(A1,1X,F11.5,1X,F11.5)
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE INTLINTOPO
+C-----------------------------------------------------------------------
+C Interpolation linéaire de la geometrie en abscisse-cote dans geomac-i0
+C Si des informations sur les compartiments sédimentaires existent,
+C on les interpole
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX,NSMAX,NCMAX,CSMAX
+ PARAMETER (LMAX=3000,NCMAX=300,NSMAX=3000,CSMAX=10)
+ INTEGER I,J,LM,II
+ INTEGER NUMCOU,NCOUAV,NCOUAP,NUMSEC,NUMCS,NBSEC
+ INTEGER NCOU(0:NSMAX),NC2(0:LMAX)
+ DOUBLE PRECISION TMAIL(LMAX),XTMAIL(LMAX),XSECT(NSMAX)
+ INTEGER NBCS(NCMAX*NSMAX)
+ DOUBLE PRECISION YCOU(NCMAX*NSMAX),ZCOU(NCMAX*NSMAX)
+ CHARACTER*1 CCOU(NCMAX*NSMAX)
+ DOUBLE PRECISION ZCS(NCMAX*NSMAX,CSMAX),DCS(NCMAX*NSMAX,CSMAX)
+ & ,SCS(NCMAX*NSMAX,CSMAX),TMCS(NCMAX*NSMAX,CSMAX)
+ INTEGER XNBCS(NCMAX*LMAX)
+ DOUBLE PRECISION XYCOU(NCMAX*LMAX),XZCOU(NCMAX*LMAX)
+ CHARACTER*1 XCCOU(NCMAX*LMAX)
+ DOUBLE PRECISION XZCS(NCMAX*LMAX,CSMAX),XDCS(NCMAX*LMAX,CSMAX)
+ & ,XSCS(NCMAX*LMAX,CSMAX),XTMCS(NCMAX*LMAX,CSMAX)
+ CHARACTER ETUDE*20,NOMFIC*40,LIGNE*630,FINLIGNE*600
+ CHARACTER*1 DESSEC
+ LOGICAL EXISTPB
+
+ COMMON/ETUDE/ETUDE
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+C Vérification: lecture du fichier mail
+C-----------------------------------------------------------------------
+ NOMFIC='mail.'//ETUDE
+ OPEN(20,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ READ(20,*)LM
+ READ(20,*)(TMAIL(I),I=1,LM)
+ READ(20,*)(XTMAIL(I),I=1,LM-1)
+ CLOSE(20)
+
+C Preambule
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'INTERPOLATION LINEAIRE DE LA TOPOGRAPHIE EN '
+ & ,'ABSCISSE-COTE'
+ WRITE(*,*)' '
+ WRITE(*,*)'ATTENTION: les sections seront interpolees point a '
+ & ,'point; il faut donc au moins qu''elles aient le meme nombre '
+ & ,'de points.'
+ WRITE(*,*)' '
+
+C Lecture de geomac-i00
+C-----------------------------------------------------------------------
+ NOMFIC='geomac-i00.'//ETUDE
+ OPEN(24,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ NCOU(0)=0
+ DESSEC=' '
+ DO NUMSEC=1,NSMAX
+C Lecture de la ligne d'en-tête de section
+ READ(24,'(A630)',END=221) LIGNE
+ READ(LIGNE,'(I4,1X,F11.3,1X,I2,1X,A600)') II,XSECT(NUMSEC)
+ & ,NCOU(NUMSEC),FINLIGNE
+cliq & ,NCOU(NUMSEC)
+ NCOU(NUMSEC)=NCOU(NUMSEC)+NCOU(NUMSEC-1)
+ IF(NUMSEC.NE.II)THEN
+ WRITE(*,*)'ATTENTION, dans ',NOMFIC,', le numero de section '
+ & ,II,' incoherent a la ',NUMSEC,'ieme maille est renumerote.'
+ ENDIF
+ IF(FINLIGNE(1:10).NE.' ')THEN
+ IF(DESSEC.EQ.' '.OR.DESSEC.EQ.'1')THEN
+ DESSEC='1'
+ ELSE
+ WRITE(*,*)'PROBLEME1 DE FORMAT SEDIMENTAIRE SECTION ',NUMSEC
+ STOP
+ ENDIF
+ ENDIF
+cliq dessec=' '
+C Lecture des paramètres sédimentaires en format simplifié
+ IF(DESSEC.EQ.'1')THEN
+ READ(FINLIGNE,*,END=222)(ZCS(NUMSEC,NUMCS),DCS(NUMSEC,NUMCS)
+ & ,SCS(NUMSEC,NUMCS),TMCS(NUMSEC,NUMCS),NUMCS=1,CSMAX)
+ 222 CONTINUE
+ NBCS(NUMSEC)=NUMCS-1
+ ENDIF
+
+ DO NUMCOU=NCOU(NUMSEC-1)+1,NCOU(NUMSEC)
+ READ(24,'(A630)') LIGNE
+ READ(LIGNE,'(A1,1X,F11.5,1X,A180)') CCOU(NUMCOU)
+ & ,YCOU(NUMCOU),FINLIGNE
+C Détermination / vérification du format
+ IF(FINLIGNE(12:22).NE.' ')THEN
+ IF(DESSEC.EQ.' '.OR.DESSEC.EQ.'2')THEN
+ DESSEC='2'
+ ELSE
+ WRITE(*,*)'PROBLEME2 DE FORMAT SEDIMENTAIRE SECTION ',NUMSEC
+ STOP
+ ENDIF
+ ELSE
+ IF(DESSEC.EQ.' '.OR.DESSEC.EQ.'0')THEN
+ DESSEC='0'
+ ELSEIF(DESSEC.EQ.'2')THEN
+ WRITE(*,*)'PROBLEME3 DE FORMAT SEDIMENTAIRE SECTION ',NUMSEC
+ STOP
+ ENDIF
+ ENDIF
+
+ IF(DESSEC.EQ.'2')THEN
+C Lecture de la ligne en format sédimentaire point par point
+ READ(FINLIGNE,*,END=224)(ZCS(NUMCOU,NUMCS),DCS(NUMCOU,NUMCS)
+ & ,SCS(NUMCOU,NUMCS),TMCS(NUMCOU,NUMCS),NUMCS=1,CSMAX)
+ 224 CONTINUE
+ NBCS(NUMCOU)=NUMCS-1
+ ELSE
+ READ(FINLIGNE,*) ZCOU(NUMCOU)
+ ENDIF
+
+ ENDDO
+ ENDDO
+ 221 CLOSE(24)
+ NBSEC=NUMSEC-1
+
+C Confirmations
+ WRITE(*,*)'Lecture des ',NBSEC,' sections de donnees de'
+ & ,' geomac-i0.',ETUDE
+ IF(DESSEC.EQ.'1')THEN
+ WRITE(*,*)'Il existe une description sedimentaire globale'
+ & ,' par section; elle sera interpolee'
+ WRITE(*,*)' '
+ ELSEIF(DESSEC.EQ.'2')THEN
+ WRITE(*,*)'Il existe une description sedimentaire'
+ & ,' point par point; elle sera interpolee'
+ WRITE(*,*)' '
+ ELSE
+ WRITE(*,*)' '
+ ENDIF
+
+C Interpolation des sections
+C-----------------------------------------------------------------------
+ WRITE(*,*)'Interpolation des sections en abscisse-cote...'
+ NC2(0)=0
+ NUMSEC=1
+ EXISTPB=.FALSE.
+ DO I=1,LM-1
+
+C Recherche des sections à interpoler
+ IF(XSECT(NUMSEC).LT.XTMAIL(I).AND.NUMSEC.LT.NBSEC)THEN
+ NUMSEC=NUMSEC+1
+ ENDIF
+
+C Calcul des couples abscisse-cote
+ IF(XSECT(NUMSEC).EQ.XTMAIL(I))THEN
+C WRITE(*,*)'Tombe pile pour I=',I
+C Retrait de ce bloc le 27jul00
+C IF(DESSEC.EQ.'1')THEN
+C XNBCS(I)=NBCS(NUMSEC)
+C DO NUMCS=1,XNBCS(I)
+C XZCS(I,NUMCS)=ZCS(NUMSEC,NUMCS)
+C XDCS(I,NUMCS)=DCS(NUMSEC,NUMCS)
+C XSCS(I,NUMCS)=SCS(NUMSEC,NUMCS)
+C XTMCS(I,NUMCS)=TMCS(NUMSEC,NUMCS)
+C ENDDO
+C ENDIF
+C Fin du bloc retiré le 27jul00
+
+ NC2(I)=NC2(I-1)+NCOU(NUMSEC)-NCOU(NUMSEC-1)
+ DO J=1,NCOU(NUMSEC)-NCOU(NUMSEC-1)
+ XYCOU(J+NC2(I-1))=YCOU(J+NCOU(NUMSEC-1))
+ XZCOU(J+NC2(I-1))=ZCOU(J+NCOU(NUMSEC-1))
+ XCCOU(J+NC2(I-1))=CCOU(J+NCOU(NUMSEC-1))
+C Ajout de ce bloc le 27jul00: description sédimentaire reconstituée
+C en chaque point si une seule donnée pour la section
+ IF(DESSEC.EQ.'1')THEN
+ XNBCS(J+NC2(I-1))=NBCS(NUMSEC)
+ DO NUMCS=1,XNBCS(J+NC2(I-1))
+ IF(NUMCS.EQ.1)THEN
+ XZCS(J+NC2(I-1),NUMCS)=XZCOU(J+NC2(I-1))
+ ELSE
+ XZCS(J+NC2(I-1),NUMCS)=XZCS(J+NC2(I-1),NUMCS-1)
+ & -(ZCS(NUMSEC,NUMCS-1)-ZCS(NUMSEC,NUMCS))
+ ENDIF
+ XDCS(J+NC2(I-1),NUMCS)=DCS(NUMSEC,NUMCS)
+ XSCS(J+NC2(I-1),NUMCS)=SCS(NUMSEC,NUMCS)
+ IF(TMCS(NUMSEC,NUMCS).LE.0.)THEN
+ XTMCS(J+NC2(I-1),NUMCS)=0.000D+00
+ ELSEIF(TMCS(NUMSEC,NUMCS).GE.999.)THEN
+ XTMCS(J+NC2(I-1),NUMCS)=999.999D+00
+ ELSE
+ XTMCS(J+NC2(I-1),NUMCS)=TMCS(NUMSEC,NUMCS)
+ ENDIF
+ ENDDO
+ ENDIF
+C Fin du bloc ajouté le 27jul00
+ IF(DESSEC.EQ.'2')THEN
+ NUMCOU=J+NCOU(NUMSEC-1)
+ XNBCS(J+NC2(I-1))=NBCS(NUMCOU)
+ DO NUMCS=1,XNBCS(J+NC2(I-1))
+ XZCS(J+NC2(I-1),NUMCS)=ZCS(NUMCOU,NUMCS)
+ XDCS(J+NC2(I-1),NUMCS)=DCS(NUMCOU,NUMCS)
+ XSCS(J+NC2(I-1),NUMCS)=SCS(NUMCOU,NUMCS)
+ IF(TMCS(NUMCOU,NUMCS).LE.0.)THEN
+ XTMCS(J+NC2(I-1),NUMCS)=0.000D+00
+ ELSEIF(TMCS(NUMCOU,NUMCS).GE.999.)THEN
+ XTMCS(J+NC2(I-1),NUMCS)=999.999D+00
+ ELSE
+ XTMCS(J+NC2(I-1),NUMCS)=TMCS(NUMCOU,NUMCS)
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDDO
+
+ ELSE
+ IF(NUMSEC.EQ.1) NUMSEC=2
+ IF((NCOU(NUMSEC)-NCOU(NUMSEC-1))
+ & .NE.(NCOU(NUMSEC-1)-NCOU(NUMSEC-2)))THEN
+ WRITE(*,*)'INTEPOLATION IMPOSSIBLE, PAS LE MÊME NOMBRE DE'
+ & ,' COUPLES ENTRE LES SECTIONS ',NUMSEC-1,' ET ',NUMSEC
+ & ,' DANS geomac-i0'
+ STOP
+ ENDIF
+C WRITE(*,*)'Interpolation pour I=',I
+C Retrait de ce bloc le 27jul00
+C IF(DESSEC.EQ.'1')THEN
+C XNBCS(I)=MIN(NBCS(NUMSEC-1),NBCS(NUMSEC))
+C DO NUMCS=1,XNBCS(I)
+C XZCS(I,NUMCS)=ZCS(NUMSEC-1,NUMCS)
+C & +(ZCS(NUMSEC,NUMCS)-ZCS(NUMSEC-1,NUMCS))
+C & *(XTMAIL(I)-XSECT(NUMSEC-1))
+C & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+C XDCS(I,NUMCS)=DCS(NUMSEC-1,NUMCS)
+C & +(DCS(NUMSEC,NUMCS)-DCS(NUMSEC-1,NUMCS))
+C & *(XTMAIL(I)-XSECT(NUMSEC-1))
+C & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+C XSCS(I,NUMCS)=SCS(NUMSEC-1,NUMCS)
+C & +(SCS(NUMSEC,NUMCS)-SCS(NUMSEC-1,NUMCS))
+C & *(XTMAIL(I)-XSECT(NUMSEC-1))
+C & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+C IF((TMCS(NUMSEC-1,NUMCS).NE.0.)
+C & .AND.(TMCS(NUMSEC,NUMCS).NE.0.))THEN
+C XTMCS(I,NUMCS)=TMCS(NUMSEC-1,NUMCS)
+C & +(TMCS(NUMSEC,NUMCS)-TMCS(NUMSEC-1,NUMCS))
+C & *(XTMAIL(I)-XSECT(NUMSEC-1))
+C & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+C ELSE
+C XTMCS(I,NUMCS)=0.
+C ENDIF
+C ENDDO
+C ENDIF
+C Fin du bloc retiré le 27jul00
+
+ NC2(I)=NC2(I-1)+NCOU(NUMSEC)-NCOU(NUMSEC-1)
+ DO J=1,NCOU(NUMSEC)-NCOU(NUMSEC-1)
+ XYCOU(J+NC2(I-1))=YCOU(J+NCOU(NUMSEC-2))
+ & +(YCOU(J+NCOU(NUMSEC-1))-YCOU(J+NCOU(NUMSEC-2)))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+ XZCOU(J+NC2(I-1))=ZCOU(J+NCOU(NUMSEC-2))
+ & +(ZCOU(J+NCOU(NUMSEC-1))-ZCOU(J+NCOU(NUMSEC-2)))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+ IF(CCOU(J+NCOU(NUMSEC-2)).EQ.CCOU(J+NCOU(NUMSEC-1)))THEN
+ XCCOU(J+NC2(I-1))=CCOU(J+NCOU(NUMSEC-2))
+ ELSE
+ WRITE(*,*)'PROBLEME DE LIMITE MINEUR/MOYEN SECTION ',I
+ XCCOU(J+NC2(I-1))='?'
+ EXISTPB=.TRUE.
+ ENDIF
+
+C Ajout de ce bloc le 27jul00: description sédimentaire reconstituée
+C en chaque point si une seule donnée pour la section
+ IF(DESSEC.EQ.'1')THEN
+ NCOUAV=NUMSEC-1
+ NCOUAP=NUMSEC
+ XNBCS(J+NC2(I-1))=MIN(NBCS(NCOUAV),NBCS(NCOUAP))
+ DO NUMCS=1,XNBCS(J+NC2(I-1))
+ IF(NUMCS.EQ.1)THEN
+ XZCS(J+NC2(I-1),NUMCS)=XZCOU(J+NC2(I-1))
+ ELSE
+ XZCS(J+NC2(I-1),NUMCS)=XZCS(J+NC2(I-1),NUMCS-1)
+ & -((ZCS(NCOUAV,NUMCS-1)-ZCS(NCOUAV,NUMCS))
+ & +((ZCS(NCOUAP,NUMCS-1)-ZCS(NCOUAP,NUMCS))
+ & -(ZCS(NCOUAV,NUMCS-1)-ZCS(NCOUAV,NUMCS)))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1)))
+ ENDIF
+ XDCS(J+NC2(I-1),NUMCS)=DCS(NCOUAV,NUMCS)
+ & +(DCS(NCOUAP,NUMCS)-DCS(NCOUAV,NUMCS))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+ XSCS(J+NC2(I-1),NUMCS)=SCS(NCOUAV,NUMCS)
+ & +(SCS(NCOUAP,NUMCS)-SCS(NCOUAV,NUMCS))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+ IF((TMCS(NCOUAV,NUMCS).GT.999.)
+ & .AND.(TMCS(NCOUAP,NUMCS).GT.999.))THEN
+ XTMCS(J+NC2(I-1),NUMCS)=999.999D+00
+ ELSEIF((TMCS(NCOUAV,NUMCS).GT.999.)
+ & .OR.(TMCS(NCOUAP,NUMCS).GT.999.))THEN
+ XTMCS(J+NC2(I-1),NUMCS)=MIN(TMCS(NCOUAV,NUMCS)
+ & ,TMCS(NCOUAP,NUMCS))
+ ELSEIF((TMCS(NCOUAV,NUMCS).NE.0.)
+ & .AND.(TMCS(NCOUAP,NUMCS).NE.0.))THEN
+ XTMCS(J+NC2(I-1),NUMCS)=TMCS(NCOUAV,NUMCS)
+ & +(TMCS(NCOUAP,NUMCS)-TMCS(NCOUAV,NUMCS))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+ ELSE
+ XTMCS(J+NC2(I-1),NUMCS)=0.000D+00
+ ENDIF
+ ENDDO
+ ENDIF
+C Fin du bloc ajouté le 27jul00
+ IF(DESSEC.EQ.'2')THEN
+ NCOUAV=J+NCOU(NUMSEC-2)
+ NCOUAP=J+NCOU(NUMSEC-1)
+ XNBCS(J+NC2(I-1))=MIN(NBCS(NCOUAV),NBCS(NCOUAP))
+ DO NUMCS=1,XNBCS(J+NC2(I-1))
+ XZCS(J+NC2(I-1),NUMCS)=ZCS(NCOUAV,NUMCS)
+ & +(ZCS(NCOUAP,NUMCS)-ZCS(NCOUAV,NUMCS))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+ XDCS(J+NC2(I-1),NUMCS)=DCS(NCOUAV,NUMCS)
+ & +(DCS(NCOUAP,NUMCS)-DCS(NCOUAV,NUMCS))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+ XSCS(J+NC2(I-1),NUMCS)=SCS(NCOUAV,NUMCS)
+ & +(SCS(NCOUAP,NUMCS)-SCS(NCOUAV,NUMCS))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+ IF((TMCS(NCOUAV,NUMCS).GT.999.)
+ & .AND.(TMCS(NCOUAP,NUMCS).GT.999.))THEN
+ XTMCS(J+NC2(I-1),NUMCS)=999.999D+00
+ ELSEIF((TMCS(NCOUAV,NUMCS).GT.999.)
+ & .OR.(TMCS(NCOUAP,NUMCS).GT.999.))THEN
+ XTMCS(J+NC2(I-1),NUMCS)=MIN(TMCS(NCOUAV,NUMCS)
+ & ,TMCS(NCOUAP,NUMCS))
+ ELSEIF((TMCS(NCOUAV,NUMCS).NE.0.)
+ & .AND.(TMCS(NCOUAP,NUMCS).NE.0.))THEN
+ XTMCS(J+NC2(I-1),NUMCS)=TMCS(NCOUAV,NUMCS)
+ & +(TMCS(NCOUAP,NUMCS)-TMCS(NCOUAV,NUMCS))
+ & *(XTMAIL(I)-XSECT(NUMSEC-1))
+ & /(XSECT(NUMSEC)-XSECT(NUMSEC-1))
+ ELSE
+ XTMCS(J+NC2(I-1),NUMCS)=0.000D+00
+ ENDIF
+ ENDDO
+ ENDIF
+
+ ENDDO
+ ENDIF
+
+ ENDDO
+
+C Ecriture de geomac-i0
+C-----------------------------------------------------------------------
+C Retrait de ce bloc le 27jul00
+C IF(DESSEC.EQ.'1')THEN
+CC Description sedimentaire globale par section
+C WRITE(24,'(I4,1X,A1)')LM-1,DESSEC
+C DO I=1,LM-1
+C WRITE(FINLIGNE,9995)(XZCS(I,NUMCS),XDCS(I,NUMCS)
+C & ,XSCS(I,NUMCS),XTMCS(I,NUMCS),NUMCS=1,XNBCS(I))
+C WRITE(24,9991)I,XTMAIL(I),NC2(I)-NC2(I-1),FINLIGNE
+C DO J=NC2(I-1)+1,NC2(I)
+C WRITE(24,9992)XCCOU(J),XYCOU(J),XZCOU(J)
+C ENDDO
+C ENDDO
+C ELSEIF(DESSEC.EQ.'2')THEN
+C Fin du bloc retiré le 27jul00 (ajout de la ligne suivante)
+ IF((DESSEC.EQ.'1').OR.(DESSEC.EQ.'2'))THEN
+ NOMFIC='geomac-i.'//ETUDE
+ WRITE(*,*)' '
+ WRITE(*,*)'Ecriture de ',NOMFIC
+ OPEN(24,FILE=NOMFIC,STATUS='UNKNOWN',FORM='FORMATTED')
+C Description sedimentaire point par point
+C WRITE(24,'(I4,1X,A1)')LM-1,DESSEC
+ WRITE(24,'(I4)')LM-1
+ DO I=1,LM-1
+ WRITE(24,9993)I,XTMAIL(I),NC2(I)-NC2(I-1)
+ DO J=NC2(I-1)+1,NC2(I)
+ WRITE(FINLIGNE,9995)(XZCS(J,NUMCS),XDCS(J,NUMCS)
+ & ,XSCS(J,NUMCS),XTMCS(J,NUMCS),NUMCS=1,XNBCS(J))
+C WRITE(24,999x)XCCOU(J),XYCOU(J),XZCOU(J),FINLIGNE
+ WRITE(24,9994)XCCOU(J),XYCOU(J),FINLIGNE
+ ENDDO
+ ENDDO
+ ELSE
+ NOMFIC='geomac-i0.'//ETUDE
+cliq NOMFIC='geomac-i.'//ETUDE
+ WRITE(*,*)' '
+ WRITE(*,*)'Ecriture de ',NOMFIC
+ OPEN(24,FILE=NOMFIC,STATUS='UNKNOWN',FORM='FORMATTED')
+C Pas de description sedimentaire
+ WRITE(24,'(I4,1X,A1)')LM-1,DESSEC
+
+ DO I=1,LM-1
+ WRITE(24,9993)I,XTMAIL(I),NC2(I)-NC2(I-1)
+ DO J=NC2(I-1)+1,NC2(I)
+ WRITE(24,9992)XCCOU(J),XYCOU(J),XZCOU(J)
+cliq WRITE(24,9996)XCCOU(J),XYCOU(J),XZCOU(J),0.1,1.,0.
+ ENDDO
+ ENDDO
+ ENDIF
+ CLOSE(24)
+
+C Traitement des erreurs
+C-----------------------------------------------------------------------
+ IF(EXISTPB)THEN
+ WRITE(*,*)' '
+ WRITE(*,*)'PROBLEME RENCONTRE LORS DE L''INTERPOLATION:'
+ WRITE(*,*)'Vous devez editer et modifier le fichier geomac-i'
+ WRITE(*,*)'Le caractère ''?'' indique l''endroit'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+ ENDIF
+ RETURN
+
+ 1000 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'ERREUR sur le fichier ',NOMFIC
+ WRITE(*,*)'Le fichier n''existe pas ou est illisible'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+ RETURN
+
+ 9991 FORMAT(I4,1X,F11.3,1X,I2,1X,A)
+ 9992 FORMAT(A1,1X,F11.5,F13.5)
+ 9993 FORMAT(I4,1X,F11.3,1X,I2)
+ 9994 FORMAT(A1,1X,F11.5,A)
+ 9995 FORMAT(10(F13.5,2F15.10,F15.5))
+ 9996 FORMAT(A1,1X,F11.5,F13.5,2F15.10,F15.5)
+
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE TRMOCRUBE
+C-----------------------------------------------------------------------
+C À partir d'un fichier de topographie au format Mocahy,
+C écrit un fichier de topographie geomac-i0,
+C crée un fichier planimétrique 'm.etude'
+C appelle l'extracteur de maillage
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX,NCMAX,NPMAX
+ PARAMETER (LMAX=3000,NCMAX=300,NPMAX=500)
+ INTEGER LM,I,J,NSECT,NPT,INSECT,INPT
+ INTEGER JAG,JAD,NTAG,NTAGPL,NTAGCG,NTAGCD,NTAGAG,NTAGAD
+ INTEGER I1,I2,I3,I4,NBTAG,NBSECT,NBPT(LMAX),INBPT(LMAX)
+ LOGICAL PTREG,PROJ,ORDPL,ORDPT,ETIG,ETID
+ CHARACTER ETUDE*20,NOMFIC*40,NOMFIC2*40,LIGNE*50,REP*1
+ CHARACTER NOMSECT(LMAX)*12,INOMSECT(LMAX)*12
+ CHARACTER TAGP(LMAX,NPMAX)*3,ITAGP(LMAX,NPMAX)*3,TAG(0:NPMAX)*3
+ CHARACTER CCOU(LMAX,NCMAX)*1
+ DOUBLE PRECISION EPS,EPSX
+ DOUBLE PRECISION XP(LMAX,NPMAX),YP(LMAX,NPMAX),ZP(LMAX,NPMAX)
+ DOUBLE PRECISION IXP(LMAX,NPMAX),IYP(LMAX,NPMAX),IZP(LMAX,NPMAX)
+ DOUBLE PRECISION XSECT(LMAX),IXSECT(LMAX),YJ(LMAX,NCMAX)
+ DOUBLE PRECISION TMAIL(LMAX),XTMAIL(LMAX)
+ DOUBLE PRECISION XINIT,XPPREC,YPPREC
+ DOUBLE PRECISION NORM,XVU,YVU
+ DOUBLE PRECISION XV(NPMAX),YV(NPMAX)
+
+ COMMON/ETUDE/ETUDE
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+C Initialisation
+C-----------------------------------------------------------------------
+C EPS est la tolérance en dessous de laquelle deux points X,Y,Z
+C sont confondus
+C EPSX est la tolérance le long du profil en long en dessous de
+C laquelle deux sections sont confondues
+ EPS=0.0001
+ EPSX=0.0001
+ TAG(0)=' '
+ WRITE(*,*)' '
+ WRITE(*,*)'SAISIE DE LA TOPOGRAPHIE A PARTIR D''UN FICHIER MOCAHY'
+ WRITE(*,*)' '
+
+C Lecture du fichier de géométrie type Mocahy *.m
+C-----------------------------------------------------------------------
+ WRITE(*,*)'Quel est le nom du fichier de depart ? (*.m)'
+ READ(*,*)NOMFIC2
+ OPEN(25,FILE=NOMFIC2,STATUS='OLD',ERR=1000)
+
+C Lecture de l'en-tête du fichier
+ 251 CONTINUE
+ READ(25,100)LIGNE
+ IF(LIGNE(1:1).EQ.'#') GOTO 251
+
+ NSECT=0
+ NBTAG=0
+ DO WHILE(.TRUE.)
+C Lecture d'un en-tête de section
+ NSECT=NSECT+1
+ NPT=0
+ READ(LIGNE,101)I1,I2,I3,I4,XSECT(NSECT),NOMSECT(NSECT)
+ PTREG=.TRUE.
+ NTAG=1
+ DO WHILE(PTREG)
+C Lecture d'un point de données
+ READ(25,100,END=252)LIGNE
+ NPT=NPT+1
+ READ(LIGNE,102)XP(NSECT,NPT),YP(NSECT,NPT),ZP(NSECT,NPT)
+ & ,TAGP(NSECT,NPT)
+C Test du dernier point de la section
+ IF((ABS(XP(NSECT,NPT)-999.9990).LE.EPS)
+ & .AND.(ABS(YP(NSECT,NPT)-999.9990).LE.EPS))THEN
+ NBPT(NSECT)=NPT-1
+ PTREG=.FALSE.
+ ENDIF
+C Test des étiquettes de lignes directrices
+ IF(TAGP(NSECT,NPT).NE.' ')THEN
+ IF(I1.EQ.1)THEN
+ TAG(NTAG)=TAGP(NSECT,NPT)
+ NBTAG=NTAG
+ NTAG=NTAG+1
+ ELSE
+ IF(TAGP(NSECT,NPT).NE.TAG(NTAG))THEN
+ WRITE(*,*)'PROBLEME avec la ligne directrice '
+ & ,TAGP(NSECT,NPT),' a la section ',NSECT
+ ENDIF
+ IF(NTAG.GT.NBTAG)THEN
+ WRITE(*,*)'PROBLEME trop de lignes directrices'
+ & ,' a la section ',NSECT
+ ENDIF
+ NTAG=NTAG+1
+ ENDIF
+ ENDIF
+ ENDDO
+ READ(25,100,END=252)LIGNE
+ ENDDO
+ 252 CONTINUE
+ NBSECT=NSECT
+C au cas ou le fichier *.m se terminerait par des lignes blanches
+ IF(NPT.EQ.0)NBSECT=NBSECT-1
+ CLOSE(25)
+
+C Affichages d'information et de contrôle sur le fichier lu
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'Lecture des ',NBSECT,' sections du fichier ',NOMFIC
+ IF(XSECT(1).LT.XSECT(NBSECT))THEN
+ ORDPL=.TRUE.
+ WRITE(*,'(A,F12.4,A,A)')' - section amont a ',XSECT(1)
+ & ,' m ',NOMSECT(1)
+ WRITE(*,'(A,F12.4,A,A)')' - section aval a ',XSECT(NBSECT)
+ & ,' m ',NOMSECT(NBSECT)
+ ELSE
+ ORDPL=.FALSE.
+ WRITE(*,'(A,F12.4,A,A)')' - section amont a ',XSECT(NBSECT)
+ & ,' m ',NOMSECT(NBSECT)
+ WRITE(*,'(A,F12.4,A,A)')' - section aval a ',XSECT(1)
+ & ,' m ',NOMSECT(1)
+ ENDIF
+ WRITE(*,*)' '
+ WRITE(*,*)'Lecture des ',NBTAG,' lignes directrices'
+ WRITE(*,'(20(1X,I3))')(NTAG,NTAG=1,NBTAG)
+ WRITE(*,'(20(1X,A3))')(TAG(NTAG),NTAG=1,NBTAG)
+
+C Récupération des informations nécessaires; profil en long
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'Les abscisses du profil en long doivent etre '
+ & ,'croissantes de l''amont vers l''aval.'
+ WRITE(*,*)'Voulez-vous changer les abscisses actuelles ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+ WRITE(*,*)'Vous avez plusieurs possibilites pour la '
+ & ,'modification des abscisses:'
+ WRITE(*,*)'- inverser l''ordre par changement de signe'
+ WRITE(*,*)'- recalculer les abscisses le long d''une directrice'
+ WRITE(*,*)'- saisir de nouvelles abscisses'
+ WRITE(*,*)' '
+ WRITE(*,*)'Voulez-vous inverser le sens en changeant le signe '
+ & ,'des abscisses en long ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+C Inversion des abscisses
+ ORDPL=.NOT.ORDPL
+ DO NSECT=1,NBSECT
+ XSECT(NSECT)=-XSECT(NSECT)
+ ENDDO
+ ELSE
+ WRITE(*,*)'Voulez-vous recalculer les abscisses en long le '
+ & ,'long d''une ligne directrice ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+C Recalcul des abscisses
+ 261 CONTINUE
+ WRITE(*,*)'Quel est le numero de cette ligne directrice ?'
+ READ(*,*)NTAGPL
+ IF((NTAGPL.LT.1).OR.(NTAGPL.GT.NBTAG))THEN
+ WRITE(*,*)'PROBLEME de ligne directrice. Recommencez.'
+ GOTO 261
+ ENDIF
+ WRITE(*,*)'Voulez-vous inverser le sens de parcours du '
+ & ,'profil en long ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+ ORDPL=.NOT.ORDPL
+ ENDIF
+ WRITE(*,*)'Quelle est l''abscisse de depart ? (m)'
+ READ(*,*)XINIT
+ IF(ORDPL)THEN
+ XSECT(1)=XINIT
+ NPT=1
+ DO WHILE(TAGP(1,NPT).NE.TAG(NTAGPL))
+ NPT=NPT+1
+ ENDDO
+ XPPREC=XP(1,NPT)
+ YPPREC=YP(1,NPT)
+ DO NSECT=2,NBSECT
+ NPT=1
+ DO WHILE(TAGP(NSECT,NPT).NE.TAG(NTAGPL))
+ NPT=NPT+1
+ ENDDO
+ XSECT(NSECT)=XSECT(NSECT-1)+((XP(NSECT,NPT)-XPPREC)**2.
+ & +(YP(NSECT,NPT)-YPPREC)**2.)**0.5
+ XPPREC=XP(NSECT,NPT)
+ YPPREC=YP(NSECT,NPT)
+ ENDDO
+ ELSE
+ XSECT(NBSECT)=XINIT
+ NPT=1
+ DO WHILE(TAGP(NBSECT,NPT).NE.TAG(NTAGPL))
+ NPT=NPT+1
+ ENDDO
+ XPPREC=XP(NBSECT,NPT)
+ YPPREC=YP(NBSECT,NPT)
+ DO NSECT=NBSECT-1,1,-1
+ NPT=1
+ DO WHILE(TAGP(NSECT,NPT).NE.TAG(NTAGPL))
+ NPT=NPT+1
+ ENDDO
+ XSECT(NSECT)=XSECT(NSECT+1)+((XP(NSECT,NPT)-XPPREC)**2.
+ & +(YP(NSECT,NPT)-YPPREC)**2.)**0.5
+ XPPREC=XP(NSECT,NPT)
+ YPPREC=YP(NSECT,NPT)
+ ENDDO
+ ENDIF
+ ELSE
+ WRITE(*,*)'Voulez-vous saisir de nouvelles abscisses en '
+ & ,'long ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+C Saisie de nouvelles abscisses
+ WRITE(*,*)'Voulez-vous inverser le sens de parcours du '
+ & ,'profil en long ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+ ORDPL=.NOT.ORDPL
+ ENDIF
+ IF(ORDPL)THEN
+ DO NSECT=1,NBSECT
+ WRITE(*,*)'Section ',NSECT,', ancienne abscisse '
+ & ,XSECT(NSECT),' m ',NOMSECT(NSECT)
+ 262 CONTINUE
+ WRITE(*,*)'Quelle est la nouvelle abscisse ? (m)'
+ READ(*,*)XSECT(NSECT)
+ IF((NSECT.GT.1)
+ & .AND.(XSECT(NSECT).LE.XSECT(NSECT-1)))THEN
+ WRITE(*,*)'ATTENTION abscisse non croissante.'
+ & ,' Recommencez.'
+ GOTO 262
+ ENDIF
+ ENDDO
+ ELSE
+ DO NSECT=NBSECT,1,-1
+ WRITE(*,*)'Section ',NBSECT-NSECT+1
+ & ,', ancienne abscisse ',XSECT(NSECT),' m '
+ & ,NOMSECT(NSECT)
+ 263 CONTINUE
+ WRITE(*,*)'Quelle est la nouvelle abscisse ? (m)'
+ READ(*,*)XSECT(NSECT)
+ IF((NSECT.LT.NBSECT)
+ & .AND.(XSECT(NSECT).LE.XSECT(NSECT+1)))THEN
+ WRITE(*,*)'ATTENTION abscisse non croissante.'
+ & ,' Recommencez.'
+ GOTO 263
+ ENDIF
+ ENDDO
+ ENDIF
+ ELSE
+ WRITE(*,*)'Bon, ben finalement, on ne fait rien...'
+ ENDIF
+ ENDIF
+ ENDIF
+ ENDIF
+
+C Remise en ordre éventuelle des sections
+C-----------------------------------------------------------------------
+ IF(.NOT.ORDPL)THEN
+C Inversion dans les variables temporaires
+ DO NSECT=1,NBSECT
+ INSECT=NBSECT-NSECT+1
+ IXSECT(INSECT)=XSECT(NSECT)
+ INOMSECT(INSECT)=NOMSECT(NSECT)
+ INBPT(NSECT)=NBPT(NSECT)
+ DO NPT=1,NBPT(NSECT)
+ IXP(INSECT,NPT)=XP(NSECT,NPT)
+ IYP(INSECT,NPT)=YP(NSECT,NPT)
+ IZP(INSECT,NPT)=ZP(NSECT,NPT)
+ ITAGP(INSECT,NPT)=TAGP(NSECT,NPT)
+ ENDDO
+ ENDDO
+
+C Retour aux variables initiales
+ DO NSECT=1,NBSECT
+ XSECT(NSECT)=IXSECT(NSECT)
+ NOMSECT(NSECT)=INOMSECT(NSECT)
+ NBPT(NSECT)=INBPT(NSECT)
+ DO NPT=1,INBPT(NSECT)
+ XP(NSECT,NPT)=IXP(NSECT,NPT)
+ YP(NSECT,NPT)=IYP(NSECT,NPT)
+ ZP(NSECT,NPT)=IZP(NSECT,NPT)
+ TAGP(NSECT,NPT)=ITAGP(NSECT,NPT)
+ ENDDO
+ ENDDO
+ WRITE(*,*)'Inversion de l''ordre des sections'
+ ENDIF
+
+C Récupération des informations nécessaires; orientation des sections
+C-----------------------------------------------------------------------
+ 270 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'Les sections en travers doivent etre parcourues de la '
+ & ,'rive gauche a la rive droite.'
+ WRITE(*,*)'Voulez-vous changer le sens des sections ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+ ORDPT=.FALSE.
+ ELSE
+ ORDPT=.TRUE.
+ ENDIF
+
+C Remise en ordre éventuelle des points des sections
+C-----------------------------------------------------------------------
+ IF(.NOT.ORDPT)THEN
+C Inversion dans les variables temporaires
+ DO NSECT=1,NBSECT
+ DO NPT=1,NBPT(NSECT)
+ INPT=NBPT(NSECT)-NPT+1
+ IXP(NSECT,INPT)=XP(NSECT,NPT)
+ IYP(NSECT,INPT)=YP(NSECT,NPT)
+ IZP(NSECT,INPT)=ZP(NSECT,NPT)
+ ITAGP(NSECT,INPT)=TAGP(NSECT,NPT)
+ ENDDO
+ ENDDO
+
+C Retour aux variables initiales
+ DO NSECT=1,NBSECT
+ DO NPT=1,NBPT(NSECT)
+ XP(NSECT,NPT)=IXP(NSECT,NPT)
+ YP(NSECT,NPT)=IYP(NSECT,NPT)
+ ZP(NSECT,NPT)=IZP(NSECT,NPT)
+ TAGP(NSECT,NPT)=ITAGP(NSECT,NPT)
+ ENDDO
+ ENDDO
+ WRITE(*,*)'Inversion de l''ordre des points des sections'
+ ENDIF
+
+C Récupération des informations nécessaires; limites du lit mineur
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'Pour un lit compose, les lignes directrices limites '
+ & ,'entre les lits mineur et moyen doivent etre choisies. Si '
+ & ,'elles ne le sont pas, seul un lit mineur est considere.'
+ 271 CONTINUE
+ WRITE(*,*)'Quel est le numero de la limite gauche ? (0 si rien)'
+ READ(*,*)NTAGCG
+ IF((NTAGCG.LT.0).OR.(NTAGCG.GT.NBTAG))THEN
+ WRITE(*,*)'PROBLEME de ligne directrice. Recommencez.'
+ GOTO 271
+ ENDIF
+ 272 CONTINUE
+ WRITE(*,*)'Quel est le numero de la limite droite ? (0 si rien)'
+ READ(*,*)NTAGCD
+ IF((NTAGCD.LT.0).OR.(NTAGCD.GT.NBTAG))THEN
+ WRITE(*,*)'PROBLEME de ligne directrice. Recommencez.'
+ GOTO 272
+ ENDIF
+ IF((NTAGCG.NE.0).AND.(NTAGCD.NE.0).AND.(NTAGCG.EQ.NTAGCD))THEN
+ WRITE(*,*)'PROBLEME lignes directrices confondues. Recommencez.'
+ GOTO 271
+ ENDIF
+ IF((NTAGCG.NE.0).AND.(NTAGCD.NE.0)
+ & .AND.(NTAGCG.GT.NTAGCD).AND.ORDPT)THEN
+ WRITE(*,*)'PROBLEME orientation de la section incompatible avec'
+ & ,' le choix des lignes directrices. Recommencez.'
+ GOTO 270
+ ENDIF
+
+C Récupération des informations nécessaires; limites du lit actif
+C-----------------------------------------------------------------------
+ 280 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'Vous pouvez choisir les lignes directrices delimitant '
+ & ,'le lit actif. Par defaut, le lit actif est confondu avec le '
+ & ,'lit mineur.'
+ WRITE(*,*)'CONSEIL: vous pouvez aussi utiliser le menu de saisie '
+ & ,'de la geometrie qui permet de rentrer les largeurs actives '
+ & ,'par zones.'
+ 281 CONTINUE
+ WRITE(*,*)'Quel est le numero de la limite gauche ? (0 si rien)'
+ READ(*,*)NTAGAG
+ IF((NTAGAG.LT.0).OR.(NTAGAG.GT.NBTAG))THEN
+ WRITE(*,*)'PROBLEME de ligne directrice. Recommencez.'
+ GOTO 281
+ ENDIF
+ 282 CONTINUE
+ WRITE(*,*)'Quel est le numero de la limite droite ? (0 si rien)'
+ READ(*,*)NTAGAD
+ IF((NTAGAD.LT.0).OR.(NTAGAD.GT.NBTAG))THEN
+ WRITE(*,*)'PROBLEME de ligne directrice. Recommencez.'
+ GOTO 282
+ ENDIF
+ IF((NTAGAG.NE.0).AND.(NTAGAD.NE.0)
+ & .AND.(NTAGAG.GT.NTAGAD).AND.ORDPT)THEN
+ WRITE(*,*)'PROBLEME orientation de la section incompatible avec'
+ & ,' le choix des lignes directrices. Recommencez.'
+ GOTO 280
+ ENDIF
+
+C Récupération des informations nécessaires; aplanissement des sections
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'La section en travers n''est pas forcement plane; '
+ & ,'pour le calcul monodimensionnel, il est possible de la '
+ & ,'projeter sur le plan vertical passant par ses extremites ou '
+ & ,'de la deplier comme un accordeon.'
+ WRITE(*,*)'Voulez-vous projeter la section ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'O').OR.(REP.EQ.'o'))THEN
+ PROJ=.TRUE.
+ ELSE
+ PROJ=.FALSE.
+ ENDIF
+
+ IF(PROJ)THEN
+C Calcul des points des sections en travers; cas section projetée
+C-----------------------------------------------------------------------
+C Le plan de projection passe par les premier
+C et dernier points de la section.
+C L'abscisse latérale YJ, qui correspond à l'abscisse latérale projetée
+C est le produit scalaire du vecteur d'un point par rapport à l'origine
+C par le vecteur de projection unitaire.
+C XV,YV coordonnées des vecteurs entre le premier point et le point n
+C XVU,YVU coordonnées du vecteur unitaire
+C entre les premier et dernier points
+ DO NSECT=1,NBSECT
+ XV(NBPT(NSECT))=XP(NSECT,NBPT(NSECT))-XP(NSECT,1)
+ YV(NBPT(NSECT))=YP(NSECT,NBPT(NSECT))-YP(NSECT,1)
+ NORM=(XV(NBPT(NSECT))**2+YV(NBPT(NSECT))**2)**0.5
+ If(norm.lT.eps)then
+ write(*,*)'attention section ',nsect
+ write(*,*)'premiers et derniers points de memes coordonnees'
+ write(*,*)'leur distance doit etre > 0.1 mm'
+ stop
+ endif
+ XVU=XV(NBPT(NSECT))/NORM
+ YVU=YV(NBPT(NSECT))/NORM
+C WRITE(*,*)' NORM=',NORM,' XVU=',XVU,' YVU=',YVU
+ YJ(NSECT,1)=0.
+ DO NPT=2,NBPT(NSECT)-1
+ XV(NPT)=XP(NSECT,NPT)-XP(NSECT,1)
+ YV(NPT)=YP(NSECT,NPT)-YP(NSECT,1)
+ YJ(NSECT,NPT)=XV(NPT)*XVU+YV(NPT)*YVU
+C WRITE(*,*)'YJ(',NSECT,',',NPT,')=',YJ(NSECT,NPT)
+C On recalcule maintenant les nouvelles coordonnées
+C (après projection sur le plan) des points de départ
+ XP(NSECT,NPT)=XP(NSECT,1)+YJ(NSECT,NPT)*XVU
+ YP(NSECT,NPT)=YP(NSECT,1)+YJ(NSECT,NPT)*YVU
+ ENDDO
+ YJ(NSECT,NBPT(NSECT))=NORM
+ ENDDO
+C Remarque: il existe un risque d'inversion de l'ordre des points
+C après projection (YJ non partout croissants)
+C Ce cas peu probable n'est pas traité dans la version actuelle,
+C mais c'est ici qu'il faudrait faire le tri des YJ,XP,YP,ZP
+C par XJ croissants
+C pas souhaitable si section presque fermee
+
+ ELSE
+C Calcul des points des sections en travers; cas section dépliée
+C-----------------------------------------------------------------------
+ DO NSECT=1,NBSECT
+ YJ(NSECT,1)=0.
+ DO NPT=2,NBPT(NSECT)
+ YJ(NSECT,NPT)=YJ(NSECT,NPT-1)
+ & +((XP(NSECT,NPT)-XP(NSECT,NPT-1))**2.
+ & +(YP(NSECT,NPT)-YP(NSECT,NPT-1))**2.)**0.5
+ ENDDO
+ ENDDO
+C fin du if sur proj
+ ENDIF
+ WRITE(*,*)' '
+ WRITE(*,*)'Calcul des sections en travers'
+
+C Élimination des informations superflues
+C-----------------------------------------------------------------------
+C NSECT (1..NBSECT),NPT (1..NBPT(NSECT)) indices
+C des informations exhaustives
+C I (1..LM-1),J (1..NBCOU) indices des informations triées
+ I=0
+ DO NSECT=1,NBSECT
+C Élimination des sections confondues
+ IF((NSECT.GT.1)
+ & .AND.(ABS(XSECT(NSECT)-XSECT(NSECT-1)).LE.EPSX))THEN
+ WRITE(*,*)'ATTENTION: la section de donnees ',NSECT-1
+ & ,' a ete supprimee car trop proche de la suivante.'
+ ELSE
+ I=I+1
+ ENDIF
+ XSECT(I)=XSECT(NSECT)
+ NOMSECT(I)=NOMSECT(NSECT)
+ J=0
+ DO NPT=1,NBPT(NSECT)
+C Élimination des points dont les abscisses latérales sont confondues...
+C modif sairube2: on elimine que les points strictement confondus
+C (en outre eps est passe de 1 mm à 0,1 mm)
+ IF((NPT.GT.1)
+ & .AND.(ABS(YJ(NSECT,NPT)-YJ(NSECT,NPT-1)).LE.EPS)
+ : .AND.(ABS(ZP(NSECT,NPT)-ZP(NSECT,NPT-1)).LE.EPS))THEN
+ WRITE(*,*)'ATTENTION: section ',NSECT,', le point ',NPT-1
+ & ,' a ete supprime car trop proche du suivant.'
+C ... mais on garde alors les étiquettes intéressantes
+ IF((TAGP(NSECT,NPT-1).EQ.TAG(NTAGCG))
+ & .AND.(TAGP(NSECT,NPT).EQ.TAG(NTAGCD)))THEN
+ TAGP(NSECT,NPT)=' '
+ ELSEIF(TAGP(NSECT,NPT-1).EQ.TAG(NTAGCG))THEN
+ TAGP(NSECT,NPT)=TAGP(NSECT,NPT-1)
+ ELSEIF(TAGP(NSECT,NPT-1).EQ.TAG(NTAGCD))THEN
+ TAGP(NSECT,NPT)=TAGP(NSECT,NPT-1)
+ ELSEIF((TAGP(NSECT,NPT-1).EQ.TAG(NTAGAG))
+ & .AND.(TAGP(NSECT,NPT).EQ.TAG(NTAGAD)))THEN
+ TAGP(NSECT,NPT)=' '
+ ELSEIF(TAGP(NSECT,NPT-1).EQ.TAG(NTAGAG))THEN
+ TAGP(NSECT,NPT)=TAGP(NSECT,NPT-1)
+ ELSEIF(TAGP(NSECT,NPT-1).EQ.TAG(NTAGAD))THEN
+ TAGP(NSECT,NPT)=TAGP(NSECT,NPT-1)
+ ENDIF
+ ELSE
+ J=J+1
+ ENDIF
+ XP(I,J)=XP(NSECT,NPT)
+ YP(I,J)=YP(NSECT,NPT)
+ ZP(I,J)=ZP(NSECT,NPT)
+ TAGP(I,J)=TAGP(NSECT,NPT)
+ YJ(I,J)=YJ(NSECT,NPT)
+ ENDDO
+ NBPT(I)=J
+ ENDDO
+ LM=I+1
+ WRITE(*,*)' '
+ WRITE(*,*)'Tri des points superflus'
+
+C Gestion des étiquettes de lignes directrices
+C-----------------------------------------------------------------------
+ DO I=1,LM-1
+C Pour TAGP, les superflues sont éliminées
+C Pour CCOU, elles sont ajoutées si elles n'existent pas
+C ETIG,ETID marques de l'existence de G et D
+ ETIG=.FALSE.
+ ETID=.FALSE.
+ DO J=1,NBPT(I)
+ IF((NTAGCG.NE.0).AND.(TAGP(I,J).EQ.TAG(NTAGCG)))THEN
+ CCOU(I,J)='G'
+ ETIG=.TRUE.
+ ELSEIF((NTAGCD.NE.0).AND.(TAGP(I,J).EQ.TAG(NTAGCD)))THEN
+ CCOU(I,J)='D'
+ ETID=.TRUE.
+ ELSEIF((NTAGAG.NE.0).AND.(TAGP(I,J).EQ.TAG(NTAGAG)))THEN
+ CCOU(I,J)='g'
+ ELSEIF((NTAGAD.NE.0).AND.(TAGP(I,J).EQ.TAG(NTAGAD)))THEN
+ CCOU(I,J)='d'
+ ELSE
+ TAGP(I,J)=' '
+ CCOU(I,J)=' '
+ ENDIF
+ ENDDO
+ IF(.NOT.ETIG) CCOU(I,1)='G'
+ IF(.NOT.ETID) CCOU(I,NBPT(I))='D'
+ ENDDO
+
+C Écriture de 'geomac-i0.etude'
+C-----------------------------------------------------------------------
+ NOMFIC='geomac-i0.'//ETUDE
+cliq NOMFIC='geomac-i.'//ETUDE
+ OPEN(10,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(10,110)LM-1
+ DO I=1,LM-1
+ WRITE(10,111)I,XSECT(I),NBPT(I)
+ DO J=1,NBPT(I)
+ WRITE(10,112)CCOU(I,J),YJ(I,J),ZP(I,J)
+cliq WRITE(10,113)CCOU(I,J),YJ(I,J),ZP(I,J),0.1,1.,0.
+ ENDDO
+ ENDDO
+ WRITE(*,*)' '
+ WRITE(*,*)'Ecriture du fichier ',NOMFIC
+ CLOSE(10)
+
+C Écriture de 'm.etude'
+C-----------------------------------------------------------------------
+ NOMFIC='m.'//ETUDE
+ OPEN(25,FILE=NOMFIC,STATUS='UNKNOWN')
+ DO I=1,LM-1
+ WRITE(25,101)I,I,0,NBPT(I),XSECT(I),NOMSECT(I)
+ DO J=1,NBPT(I)
+ WRITE(25,102)XP(I,J),YP(I,J),ZP(I,J),TAGP(I,J)
+ ENDDO
+ WRITE(25,102)999.9990,999.9990,0.0000,' '
+ ENDDO
+ CLOSE(25)
+ WRITE(*,*)' '
+ WRITE(*,*)'Ecriture du fichier ',NOMFIC
+
+C Écriture du fichier de maillage 'mail.etude'
+C-----------------------------------------------------------------------
+ NOMFIC='geomac-i0.'//ETUDE
+cliq NOMFIC='geomac-i.'//ETUDE
+ CALL EXTMAIL(NOMFIC)
+ WRITE(*,*)' '
+ WRITE(*,*)'Ecriture du fichier mail.',ETUDE
+
+C Écriture de 'litactif.etude'
+C-----------------------------------------------------------------------
+C JAG,JAD indices des marques de lit actif à gauche 'g' et à droite 'd'
+ NOMFIC='litactif.'//ETUDE
+ OPEN(20,FILE=NOMFIC,STATUS='UNKNOWN')
+ DO I=1,LM-1
+ JAG=1
+ ETIG=.TRUE.
+ JAD=NBPT(I)
+ ETID=.TRUE.
+ DO J=1,NBPT(I)
+ IF(ETIG)THEN
+ IF(CCOU(I,J).EQ.'g')THEN
+ JAG=J
+ ETIG=.FALSE.
+ ELSEIF(CCOU(I,J).EQ.'G')THEN
+ JAG=J
+ ENDIF
+ ENDIF
+ IF(ETID)THEN
+ IF(CCOU(I,J).EQ.'d')THEN
+ JAD=J
+ ETID=.FALSE.
+ ELSEIF(CCOU(I,J).EQ.'D')THEN
+ JAD=J
+ ENDIF
+ ENDIF
+ ENDDO
+C WRITE(20,120) XSECT(I),ABS(YJ(I,JAD)-YJ(I,JAG))
+ WRITE(20,120) XSECT(I)
+ :,ABS(0.5*(YJ(I,JAD)+YJ(I,JAD-1))-0.5*(YJ(I,JAG)+YJ(I,JAG+1)))
+ ENDDO
+ CLOSE(20)
+ WRITE(*,*)' '
+ WRITE(*,*)'Ecriture du fichier ',NOMFIC
+
+ RETURN
+
+C Traitement des erreurs
+C-----------------------------------------------------------------------
+ 1000 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'ERREUR sur le fichier ',NOMFIC2
+ WRITE(*,*)'Le fichier n''existe pas ou est illisible'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+
+ 100 FORMAT(A50)
+ 101 FORMAT(4I6,F13.4,1X,A42)
+ 102 FORMAT(3f13.4,1X,A3)
+ 110 FORMAT(I4)
+ 111 FORMAT(I4,1X,F11.3,1X,I4)
+ 112 FORMAT(A1,1X,F11.5,F13.5)
+ 113 FORMAT(A1,1X,F11.5,F13.5,2F15.10,f15.5)
+ 120 FORMAT(1X,F11.3,1X,F11.5)
+
+ RETURN
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE SAISED
+C-----------------------------------------------------------------------
+C Saisie des paramètres sédimentaires par tronçon pour compléter
+C 'geomac-i0.etude' en 'geomac-i.etude'
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX,NCMAX,CSMAX
+ PARAMETER(LMAX=3000,NCMAX=300,CSMAX=10)
+ LOGICAL ENCOREZONE,ENCORECS,EPCS(LMAX,CSMAX)
+ INTEGER I,N,NBSECT,NZONE,NBZONE
+ INTEGER J,NCOU(LMAX)
+ INTEGER K,NCS,NBCS(LMAX),NBCSZONE(LMAX)
+ CHARACTER ETUDE*20,NOMFIC*40,FINLIGNE*600,REP*1
+ CHARACTER*1 CCOU(LMAX,NCMAX)
+ DOUBLE PRECISION EPSX,ZERO
+ DOUBLE PRECISION YCOU(LMAX,NCMAX),ZCOU(LMAX,NCMAX)
+ DOUBLE PRECISION XSECT(LMAX),ZCS(LMAX,NCMAX,CSMAX)
+ & ,DCS(LMAX,NCMAX,CSMAX),SCS(LMAX,NCMAX,CSMAX)
+ & ,TMCS(LMAX,NCMAX,CSMAX)
+ DOUBLE PRECISION XZONE(LMAX),ZEZONE(LMAX,CSMAX),DZONE(LMAX,CSMAX)
+ & ,SZONE(LMAX,CSMAX),TMZONE(LMAX,CSMAX)
+
+ COMMON/ETUDE/ETUDE
+
+C Initialisation
+C-----------------------------------------------------------------------
+C EPSX est la tolérance le long du profil en long en dessous de
+C laquelle deux sections sont confondues
+ EPSX=0.001
+ ZERO=0.0000000000
+ NOMFIC='geomac-i0.'//ETUDE
+ WRITE(*,*)' '
+ WRITE(*,*)'SAISIE DES PARAMETRES SEDIMENTAIRES'
+ WRITE(*,*)'POUR COMPLEMENT DU FICHIER DE TOPOGRAPHIE ',NOMFIC
+ WRITE(*,*)' '
+
+C Lecture de 'geomac-i0.etude'
+C-----------------------------------------------------------------------
+ OPEN(24,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ READ(24,*,ERR=1001)NBSECT
+ DO I=1,NBSECT
+ READ(24,*,ERR=1001)N,XSECT(I),NCOU(I)
+ IF(N.NE.I) WRITE(*,*)'ATTENTION, la section ',N
+ & ,'est mal ordonee dans ',NOMFIC
+ DO J=1,NCOU(I)
+ READ(24,130,ERR=1001)CCOU(I,J),FINLIGNE
+ READ(FINLIGNE,*,ERR=1001)YCOU(I,J),ZCOU(I,J)
+ ENDDO
+ ENDDO
+ CLOSE(24)
+C WRITE(*,*)'Lecture des ',NBSECT,' sections de donnees de ',NOMFIC
+
+C Saisie des donnee sedimentaires par zones
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'ATTENTION1: les abscisses des limites de zones doivent'
+ & ,' etre croissantes.'
+ WRITE(*,*)'ATTENTION2: les caracteristiques sedimentaires sont '
+ & ,'constantes sur une zone (pas d''interpolation).'
+ WRITE(*,*)'ATTENTION3: dans une zone, il y a au maximum ',CSMAX
+ & ,' couches sedimentaires de caracteristiques differentes.'
+ ENCOREZONE=.TRUE.
+ NZONE=1
+ DO WHILE(ENCOREZONE)
+
+C Abscisse de fin de zone
+ WRITE(*,*)' '
+ IF(NZONE.GT.1) WRITE(*,*)'On est dans la zone ',NZONE
+ & ,' commencant a l''abscisse ',XZONE(NZONE-1),' m (exclus)'
+ 100 CONTINUE
+ WRITE(*,*)'Quelle est l''abscisse de fin de la zone ',NZONE
+ & ,'(inclus; prolonge en aval si derniere zone) ? (m)'
+ READ(*,*)XZONE(NZONE)
+ IF((NZONE.GT.1).AND.(XZONE(NZONE).LE.XZONE(NZONE-1)))THEN
+ WRITE(*,*)'ATTENTION: abscisse non croissante. Recommencez.'
+ GOTO 100
+ ENDIF
+
+C Informations pour la strate superficielle
+ NCS=1
+ WRITE(*,*)' '
+ WRITE(*,*)'On est dans la strate 1 (superficielle) '
+ & ,'de la zone',NZONE,'.'
+ WRITE(*,*)'- diametre representatif (m) ?'
+ READ(*,*)DZONE(NZONE,NCS)
+ WRITE(*,*)'- parametre d''etendue (m/m) ? [racine(D84/D16)]'
+ READ(*,*)SZONE(NZONE,NCS)
+ WRITE(*,*)'- contrainte critique de mise en mouvement (N/m2) '
+ & ,'? [0.000 si calcul par Shields; 999.999 si inerodable]'
+ READ(*,*)TMZONE(NZONE,NCS)
+ IF(TMZONE(NZONE,NCS).GE.999.)THEN
+ TMZONE(NZONE,NCS)=999.9990D+00
+ ENCORECS=.FALSE.
+ ELSE
+ IF(TMZONE(NZONE,NCS).LE.ZERO) TMZONE(NZONE,NCS)=0.D+00
+ WRITE(*,*)'Dans cette zone, y a-t-il une autre strate '
+ & ,'sedimentaire plus profonde ? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ ENCORECS=.TRUE.
+ EPCS(NZONE,NCS)=.TRUE.
+ NCS=NCS+1
+ ELSE
+ ENCORECS=.FALSE.
+ ENDIF
+ ENDIF
+
+ DO WHILE(ENCORECS)
+C Type de caractérisation de la strate: epaisseur ou cote plancher
+C EPCS(NZONE,NCS) et ZEZONE(NZONE,NCS) caractérisent la strate NCS-1
+C (son épaisseur ou son plancher) mais n'ont de sens
+C que si la strate NCS existe
+C EPCS(NZONE,1) et ZEZONE(NZONE,1) n'ont donc pas
+C vraiment de sens physique
+ IF(EPCS(NZONE,NCS-1))THEN
+ WRITE(*,*)'La strate superieure peut etre caracterisee par '
+ & ,'une epaisseur constante sur le profil en travers ou '
+ & ,'par une cote de plancher (plafond de la strate '
+ & ,'inferieure) constante.'
+ WRITE(*,*)'Voulez-vous saisir une epaisseur ? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ EPCS(NZONE,NCS)=.TRUE.
+ ELSE
+ EPCS(NZONE,NCS)=.FALSE.
+ ENDIF
+ ELSE
+ EPCS(NZONE,NCS)=.FALSE.
+ ENDIF
+
+C Saisie epaisseur ou cote plancher
+ IF(EPCS(NZONE,NCS))THEN
+ WRITE(*,*)'- epaisseur de la strate sedimentaire superieure'
+ & ,' (m) ?'
+ ELSE
+ WRITE(*,*)'- cote du plancher de la strate sedimentaire '
+ & ,'superieure (m) ?'
+ ENDIF
+ READ(*,*)ZEZONE(NZONE,NCS)
+
+C Saisie de la strate sous-jacente: diametre, étendue, contrainte
+ WRITE(*,*)' '
+ WRITE(*,*)'On est dans la strate ',NCS,' de la zone',NZONE,'.'
+ WRITE(*,*)'- diametre representatif (m) ?'
+ READ(*,*)DZONE(NZONE,NCS)
+ WRITE(*,*)'- parametre d''etendue (m/m) ? [racine(D84/D16)]'
+ READ(*,*)SZONE(NZONE,NCS)
+ WRITE(*,*)'- contrainte critique de mise en mouvement (N/m2) '
+ & ,'? [0.000 si calcul par Shields; 999.999 si inerodable]'
+ READ(*,*)TMZONE(NZONE,NCS)
+
+C Passage à la saisie du compartiment sédimentaire suivant
+ IF(TMZONE(NZONE,NCS).GE.999.)THEN
+ TMZONE(NZONE,NCS)=999.9990D+00
+ ENCORECS=.FALSE.
+ ELSEIF(NCS.EQ.CSMAX)THEN
+ ENCORECS=.FALSE.
+ ELSE
+ IF(TMZONE(NZONE,NCS).LE.ZERO) TMZONE(NZONE,NCS)=0.D+00
+ WRITE(*,*)'Dans cette zone, y a-t-il une autre strate '
+ & ,'sedimentaire plus profonde ? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ ENCORECS=.TRUE.
+ NCS=NCS+1
+ ELSE
+ ENCORECS=.FALSE.
+ ENDIF
+ ENDIF
+ ENDDO
+ NBCSZONE(NZONE)=NCS
+
+C Passage à la saisie de la zone suivante
+ WRITE(*,*)' '
+ WRITE(*,*)'Y a-t-il une autre zone ? (O/N)'
+ READ(*,*)REP
+ IF(REP.EQ.'O'.OR.REP.EQ.'o')THEN
+ ENCOREZONE=.TRUE.
+ NZONE=NZONE+1
+ ELSE
+ ENCOREZONE=.FALSE.
+ ENDIF
+ ENDDO
+ NBZONE=NZONE
+
+C Interpolation des caractéristiques sédimentaires
+C-----------------------------------------------------------------------
+ NZONE=1
+ DO I=1,NBSECT
+C Recherche de la zone de données à exploiter
+ DO WHILE((NZONE.LT.NBZONE).AND.(XZONE(NZONE).LE.XSECT(I)-EPSX))
+ NZONE=NZONE+1
+ ENDDO
+C Interpolation en chaque point
+ DO J=1,NCOU(I)
+ DO K=1,NBCSZONE(NZONE)
+ IF(K.EQ.1)THEN
+ ZCS(I,J,K)=ZCOU(I,J)
+ ELSE
+ IF(EPCS(NZONE,K))THEN
+ ZCS(I,J,K)=ZCS(I,J,K-1)-ZEZONE(NZONE,K)
+ ELSE
+ ZCS(I,J,K)=ZEZONE(NZONE,K)
+ IF(ZCS(I,J,K).GT.ZCS(I,J,K-1)) ZCS(I,J,K)=ZCS(I,J,K-1)
+ ENDIF
+ ENDIF
+ DCS(I,J,K)=DZONE(NZONE,K)
+ SCS(I,J,K)=SZONE(NZONE,K)
+ TMCS(I,J,K)=TMZONE(NZONE,K)
+ ENDDO
+ ENDDO
+ NBCS(I)=NBCSZONE(NZONE)
+ ENDDO
+ WRITE(*,*)'Fin de l''interpolation'
+
+C Écriture de 'geomac-i.etude'
+C-----------------------------------------------------------------------
+ NOMFIC='geomac-i.'//ETUDE
+ OPEN(24,FILE=NOMFIC,STATUS='UNKNOWN')
+ WRITE(24,110)NBSECT,0.
+ DO I=1,NBSECT
+ WRITE(24,111)I,XSECT(I),NCOU(I)
+ DO J=1,NCOU(I)
+ DO K=1,NBCS(I)
+ N=1+58*(K-1)
+ WRITE(FINLIGNE(N:),113)ZCS(I,J,K),DCS(I,J,K),SCS(I,J,K)
+ & ,TMCS(I,J,K)
+ ENDDO
+ WRITE(24,112)CCOU(I,J),YCOU(I,J),FINLIGNE
+ ENDDO
+ ENDDO
+ CLOSE(24)
+ WRITE(*,*)'Ecriture de ',NOMFIC
+ RETURN
+
+C Traitement des erreurs
+C-----------------------------------------------------------------------
+ 1000 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'ERREUR sur le fichier ',NOMFIC
+ WRITE(*,*)'Le fichier n''existe pas ou est illisible'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+ RETURN
+
+ 1001 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'ERREUR sur le fichier ',NOMFIC
+ WRITE(*,*)'Le format du fichier est incorrect'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+ RETURN
+
+
+
+ 110 FORMAT(I4,1X,F12.3)
+ 111 FORMAT(I4,1X,F11.3,1X,I4)
+ 112 FORMAT(A1,1X,F11.5,A600)
+ 113 FORMAT(F13.5,F15.10,F15.10,F15.5)
+ 130 FORMAT(A1,1X,A)
+ END
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE SAILACT
+C-----------------------------------------------------------------------
+C Saisie du fichier des largeurs actives 'litactif.etude'
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX
+ PARAMETER (LMAX=3000)
+ INTEGER I,NBSECT
+ LOGICAL ENCORE
+ CHARACTER ETUDE*20,NOMFIC*40,REP*1
+ DOUBLE PRECISION XSECT(LMAX),LACT(LMAX)
+
+ COMMON/ETUDE/ETUDE
+
+C Initialisation
+C-----------------------------------------------------------------------
+ NOMFIC='litactif.'//ETUDE
+ ENCORE=.TRUE.
+ I=0
+ WRITE(*,*)' '
+ WRITE(*,*)'SAISIE DES LARGEURS ACTIVES PAR ZONES'
+ WRITE(*,*)' '
+
+C Test de l'existence préalable du fichier
+C-----------------------------------------------------------------------
+ OPEN(20,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ WRITE(*,*)'Le fichier de largeurs actives ',NOMFIC
+ WRITE(*,*)'existe deja. Il a pu etre cree a partir du fichier '
+ & ,'de topographie au format Mocahy.'
+ WRITE(*,*)'Voulez-vous reellement l''ecraser ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'N').OR.(REP.EQ.'n')) RETURN
+ CLOSE(20)
+ 1000 CONTINUE
+
+C Saisie des largeurs actives par zones
+C-----------------------------------------------------------------------
+ WRITE(*,*)' '
+ WRITE(*,*)'ATTENTION: Les abscisses doivent etre croissantes.'
+ DO WHILE(ENCORE)
+ I=I+1
+ WRITE(*,*)' '
+ 1001 CONTINUE
+ WRITE(*,*)'Entrez le couple abscisse (m) - largeur active (m) '
+ & ,'numero ',I
+ READ(*,*)XSECT(I),LACT(I)
+ IF((I.GT.1).AND.(XSECT(I).LE.XSECT(I-1)))THEN
+ WRITE(*,*)'ATTENTION: abscisse non croissante. Recommencez.'
+ GOTO 1001
+ ENDIF
+ IF(LACT(I).LT.0.)THEN
+ WRITE(*,*)'ATTENTION: la largeur active ne peut pas etre '
+ & ,'negative. Recommencez.'
+ GOTO 1001
+ ENDIF
+ WRITE(*,*)'Y-a-il un autre couple ? (O/N)'
+ READ(*,'(A1)')REP
+ IF((REP.EQ.'N').OR.(REP.EQ.'n')) ENCORE=.FALSE.
+ ENDDO
+ NBSECT=I
+
+C Écriture de 'litactif.etude'
+C-----------------------------------------------------------------------
+ OPEN(20,FILE=NOMFIC,STATUS='UNKNOWN')
+ DO I=1,NBSECT
+ WRITE(20,120)XSECT(I),LACT(I)
+ ENDDO
+ CLOSE(20)
+ WRITE(*,*)' '
+ WRITE(*,*)'Ecriture du fichier ',NOMFIC
+
+ RETURN
+ 120 FORMAT(1X,F11.3,1X,F11.5)
+ END
+
+
+
+C-----------------------------------------------------------------------
+ SUBROUTINE EXTMAIL(NOMFIC)
+C-----------------------------------------------------------------------
+C Extraction du maillage à partir du fichier de géométrie NOMFIC
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ INTEGER LMAX
+ PARAMETER(LMAX=3000)
+ CHARACTER NOMFIC*38
+ INTEGER LM
+ INTEGER I,J,N,NBSEC,NBPT
+ DOUBLE PRECISION TMAIL(LMAX),XTMAIL(LMAX)
+ DOUBLE PRECISION EPSX,XSECT(LMAX)
+
+ COMMON/TMAIL/TMAIL,XTMAIL
+ COMMON/NBMA/LM
+
+C Initialisation
+C-----------------------------------------------------------------------
+C EPSX est la tolérance en espace le long du profil en long
+C en dessous de laquelle deux sections sont confondues
+ EPSX=0.01
+
+C Lecture du fichier de géométrie
+C-----------------------------------------------------------------------
+ OPEN(11,FILE=NOMFIC,STATUS='OLD',ERR=1000)
+ READ(11,*)NBSEC
+ DO I=1,NBSEC
+ READ(11,*,ERR=1001) N,XSECT(I),NBPT
+ IF(N.NE.I) WRITE(*,*)'ATTENTION, section ',N
+ & ,' mal numerotee dans le fichier ',NOMFIC
+ DO J=1,NBPT
+ READ(11,*,ERR=1001)
+ ENDDO
+ ENDDO
+ CLOSE(11)
+
+C Interpolation du maillage
+C-----------------------------------------------------------------------
+ XTMAIL(1)=XSECT(1)
+ TMAIL(1)=XSECT(1)
+ DO I=2,NBSEC
+ XTMAIL(I)=XSECT(I)
+ TMAIL(I)=0.5*(XTMAIL(I)+XTMAIL(I-1))
+ ENDDO
+ LM=NBSEC+1
+ TMAIL(LM)=XTMAIL(LM-1)
+
+C Écriture du fichier mail.ETUDE
+C-----------------------------------------------------------------------
+ CALL ECRIMA
+ RETURN
+
+C Traitement des erreurs
+C-----------------------------------------------------------------------
+ 1000 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'ERREUR sur le fichier ',NOMFIC
+ WRITE(*,*)'Le fichier n''existe pas ou est illisible'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+ RETURN
+
+ 1001 CONTINUE
+ WRITE(*,*)' '
+ WRITE(*,*)'ERREUR sur le fichier ',NOMFIC
+ WRITE(*,*)'Le format du fichier est incorrect'
+ WRITE(*,*)'Tapez <entree>'
+ READ(*,*)
+ RETURN
+
+ END
+
--
GitLab