source: ether_geisa/trunk/pgm03/analy.f @ 848

Last change on this file since 848 was 1, checked in by cbipsl, 18 years ago

Geisa inital import

File size: 11.0 KB
Line 
1C     CETTE OPTION DOIT VENIR APRES LA MODIFICATION PAR
2C     &GEISA PGM='TRS',UNITE='BINAIRE',MODIF='OUI' &END
3C     LE PROGRAMME ANL MODIFIE LA BANQUE DE LA FACON SUIVANTE POUR
4C     STOCKER LES FREQUENCES DES ISOTOPES DES MOLECULES :
5C     EN RECORD 1 AJOUTER APRES LL1 LA VALEUR  DE LL2=1
6C     A PARTIR DU RECORD IFIN+LL1 ECRIRE LES VALEURS SUIVANTES :
7C     NTAB,(TAB(J),J=1,NTAB),(YMOY(J),J=1,20),(YMAX(J),J=1,20),
8C     (ALF(J),J=1,20)
9C     MODE=-1  APPEL NORMAL DE ANL POUR LISTER LES FREQUENCES MOLECULES
10C              ENTRE NU1 ET NU2
11C     MODE=0   MODIFICATION DE LA BANQUE (VOIR PRECEDEMMENT) DANS CE CAS
12C              MODIF='OUI'
13C     MODE=1   LISTE DES FREQUENCES  DE LA BANQUE PAR MOLECULE SANS
14C              LECTURE DU FICHIER(OPTION PROVENANT DE PGM='INF')
15C       
16C     ANALYSE PAR MOLECULE ET PAR VARIETE ISOTOPIQUE DU CONTENU
17C     DE LA BANQUE DANS UN DOMAINE SPECTRAL DONNE
18C     HISTOGRAMME DES FREQUENCES DES RAIES POUR UNE CLASSE
19C     D'INTENSITE (K=1) OU DE NIVEAU DE BASE (K=2)
20C       
21C     NU1,NU2 LIMITES INF ET SUP DU DOMAINE SPECTRAL ETUDIE
22C     DNU PAS D'ETUDE(JUSQU'A 500 PAS SONT PREVUS)
23C     PAR DEFAUT TOUT LE FICHIER AVEC 1 PAS D'ETUDE
24C     XMIN(K),XMAX(K),PAS(K) RESPECTIVEMENT MINIMUM,MAXIMUM,LARGEUR
25C     DE LA CLASSE POUR LES VARIABLES INTENSITE OU NIVEAU DE BASE
26C     TAB(I) TABLEAU INDICE PAR LA VALEUR DES ISOTOPES ET CONTENANT LES
27C            FREQUENCES DES ISOTOPES
28C     CODE(I) TABLEAU CONTENANT LES CODES MOLECULES
29C     NN      TABLEAU CONTENANT LE NOMBRE D'ISOTOPES PAR MOLECULE
30C             SUIVI DES CODES ISOTOPES
31C     ALF(I)  MOYENNE DES 1/2 LARGEUR A MI-HAUTEUR PAR MOLECULE
32C     HIST(I,J) TABLEAU CONTENANT L'HISTOGRAMME A REPRESENTER
33C               I=NUMERO DE LA CLASSE
34C               J=1 HISTOGRAMME DES INTENSITES
35C               J=3 HAUTEUR DES CLASSES POUR LES INTENSITES
36C               J=2 HISTOGRAMME DES NIVEAUX DE BASE
37C               J=4 HAUTEUR DES CLASSES POUR LES NIVEAUX DE BASE
38C     ICAR  HAUTEUR MAX D'UNE CLASSE(45 CARACTERES)
39C     IVAL=0 ANALYSE POUR TOUTES MOLECULES ET TOUS ISOTOPES(OPTION PAR D
40C     IVAL#0 ANALYSE POUR UNE MOLECULE OU UN ISOTOPE
41C            IVAL=NUMERO DE LA MOLECULE OU DE L'ISOTOPE
42C     IRAN=1 POUR MOLECULE
43C     IRAN=2 POUR ISOTOPE
44C     MOLE=MOLECULE DEMANDEE EXEMPLE MOLE='H2O' OU 'CO2' ...
45C     ISOT=ISOTOPE DEMANDE EXEMPLE ISOT=161 OU 162 ...
46C     ANAL='OUI' ANALYSE PAR MOLECULE
47C     ANAL=AUTRE PAS D'ANALYSE
48C     HISTO='INTE' HISTOGRAMME INTENSITE (KHIST=1)
49C     HISTO='BASE' HISTOGRAMME NIVEAU DE BASE (KHIST=2)
50C     HISTO='DEUX' HISTOGRAMME INTENSITE ET NIVEAU DE BASE (KHIST=-1)
51C     HISTO=AUTRE PAS D'HISTOGRAMME (OPTION PAR DEFAUT)
52C     NBCLAS NOMBRE DE CLASSES DE L'HISTOGRAMME (PAR DEFAUT 10)
53C     NHIST NOMBRE MAXIMUM DE CLASSES(50)
54C     KANAL=0        HISTO
55C     KANAL=1        ANAL
56C     IUNI UNITE LOGIQUE CORRESPONDANT AU FICHIER SPECTRAL
57C       
58C --> SUBROUTINES APPELEES : IMPANL
59C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *
60C  LAST MODIF:  06.05.91 PASSAGE A 75 MOLECULES DANS LES COMMON
61C       
62C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *
63      subroutine analy(tab,hist,xmin,xmax,ymoy,ymax,alf,qq,*)
64C       
65      logical*1 qq(1)
66      character*9 trs1,trs2
67      character*7 form,bin
68      character*4 code,ctlg,mole,blanc,base,deux
69      character*3 pgm,ianl,iext,itrs,ilst,icop,info,icre,
70     &            oui,liste,iopt,modif,trans
71      character*2 ikod,icod
72      character*1 bl
73      dimension pas(2),kode(2),alf(40)
74C  NOUVEAU TABLEAU V DE LECTURE ON PASSE DE 16 -> 29 OCTETS
75      real*8 ymoyd(40),ymaxd(40),cor,a2
76      real v(29),ymoy(40),ymax(40),xmin(40),xmax(40)
77      real nu1,nu2
78      integer hist(1),tab(1),vers
79C       
80      common/p1/ nu1,nu2,dnu,nbclas,khist,kanal,vers,mode,liste
81      common/p2/ mole(75),isot(100),nbi(75),nbt(75),iopt,ctlg,trans,
82     &           trs1,trs2
83      common/p3/ imole,iran,ival,pgm,ianl,iext,itrs,ilst,icop,info,icre
84      common/p4/ nmol,knmol,ksot,kksot,ntab,nhist,kp,lre,form,bin,modif
85      common/p5/ code(75),nn(150),nq(150),ikod(18),icod(18),jdeb(75)
86      common/inteh/ incr,pas1,pmax
87      common/entsor/iuni,juni,kuni,ient,isor,iper,nresv,maxx,blanc,oui
88C       
89      equivalence (kode(1),izot),(kode(2),imol)
90      equivalence (a,v(1)),(izot,v(15)),(imol,v(16))
91C       
92      data nbklas/500/,ibase/950/,cor/1.d50/
93C       
94C     INITIALISATION DES TABLEAUX
95C       
96c     print *,' entree dans ANL'
97      impr=0
98c     print *,' mode=',mode
99      if(mode.eq.-1) go to 6
100c     print *,' appel PGEISA'
101      call pgeisa(0.,99999.,999)
102c     print *,' retour PGEISA'
103999   read (iuni,rec=1)
104     &nu1,nu2,anu,n203,nbraie,nbmol,iecr,ifin,ll1,ll2,ll3,ll4
105      if(mode.eq.0) go to 7
106      if(mode.eq.1.and.ll2.ne.0) go to 5
107      write(isor,2000)
1082000  format(///' *inf*   cette option est uniquement valable pour la ba
109     &nque'/9x,'des donnees spectroscopiques'///)
110      go to 900
1115     continue
112      ifin=ifin+ll1
113      read (iuni,rec=ifin)
114     &ntab,(tab(j),j=1,ntab),nmol,(ymoy(j),j=1,nmol),
115     &(ymax(j),j=1,nmol),(alf(j),j=1,nmol)
116      do 12 j=1,nmol
117      ymoyd(j)=ymoy(j)*(1./cor)
118      ymaxd(j)=ymax(j)*(1./cor)
11912    continue
120      call impanl(tab,nu1,nu2,impr,ymoyd,ymaxd,alf,qq,ibase)
121      go to 900
1226     continue
123      call pgeisa(nu1,nu2,*900)
124      read (iuni,rec=1)
125     &rnu1,rnu2,anu,n203,nbraie,nbmol,iecr,ifin,ll1,ll2,ll3,ll4
126      if(nu1.eq.0) nu1=rnu1
127      if(nu2.eq.99999) nu2=rnu2
1287     continue
129      do 1 j=1,nmol
130      ymoy(j)=0.
131      ymax(j)=0.
132      ymoyd(j)=0.
133      ymaxd(j)=0.
1341     alf(j)=0.
135      nbc2=nbclas*2
136      npas=nbklas
137      do 10 j=1,npas
138      xmin(j)=1.e20
139      xmax(j)=0.
14010    continue
141      lk=ntab+nbc2
142      do 11 j=1,lk
14311    tab(j)=0
144C       
145C     CALCUL DES FREQUENCES PAR MOLECULE ET PAR ISOTOPE ET DES
146C     MOYENNES ALFMOY
147C       
148      if(dnu.ne.0.) im=(nu2-nu1)/dnu
149      if(im.gt.nbklas) nu2=nu1+float(nbklas)*dnu
150      im=0
151      anu1=nu1
152      anu2=nu1+dnu
153      if(dnu.eq.0.) anu2=nu2
154100   continue
155      call lgeisa(v,*200)
156C       
157C     TRAITEMENT PARTICULIER DES MOLECULES AYANT DES ISOTOPES DUPLIQUES
158C     MOLE= 9 10 17 18 19 25 26 28 30 31 36 42
159C       
160C     DANS LE TABLEAU TAB LES INDICES  951 A 992 SONT RESERVES POUR
161C     LES ISOTOPES DUPLIQUES (IBASE=950)
162C BB 06.05.97 cas c2h4 ( 2 isotopes dupliques le 2emendice=ibase=950
163C       
164C            H2O  CO2  O3  N2O  CO   CH4  O2   NO   SO2  NO2  NH3  PH3
165      go to (541, 541,541, 541, 541, 541, 541, 541, 509, 510, 541, 541,
166C            HNO3 OH   HF   HCL  HBR  HI   CLO  OCS  H2CO C2H6 CH3D C2H2
167     &       541, 541, 541, 541, 517, 518, 519, 541, 541, 541, 541, 541,
168C            C2H4 GEH4 HCN  C3H8 C2N2 C4H2 HC3N HOCL N2   CH3CLH2O2 H2S
169     &       525, 526, 541, 528, 541, 530, 531, 541, 541, 541, 541, 536,
170C            HCOOH COF2 SF6 C3H4  HO2  CLONO2
171     &       541,  541, 541, 541, 537,  541 ),imol
172C       
173C     SO2
174509   continue
175      if(izot.eq.626) izot=ibase+imol
176      go to 541
177C       
178C     NO2
179510   continue
180      if(izot.eq.646) izot=ibase+imol
181      go to 541
182C       
183C     HBR
184517   continue
185      if(izot.eq. 19) izot=ibase+imol
186      go to 541
187C       
188C     HI
189518   continue
190      if(izot.eq. 17) izot=ibase+imol
191      go to 541
192C       
193C     CLO
194519   continue
195      if(izot.eq. 56) izot=ibase+imol
196      go to 541
197C       
198C     C2H4
199525   continue
200      if(izot.eq.211) izot=ibase+imol
201      if(izot.eq.311) izot=ibase
202      go to 541
203C       
204C     GEH4
205526   continue
206      if(izot.eq.411) izot=ibase+imol
207      go to 541
208C       
209C     C3H8
210528   continue
211      if(izot.eq.221) izot=ibase+imol
212      go to 541
213C       
214C     C4H2
215530   continue
216      if(izot.eq.211) izot=ibase+imol
217      go to 541
218C       
219C     HC3N
220531   continue
221      if(izot.eq.124) izot=ibase+imol
222      go to 541
223C       
224C     H2S
225536   continue
226      if(izot.eq.131) izot=ibase+imol
227C     GO TO 541
228C
229C     HO2
230537   continue
231      if(izot.eq.166) izot=ibase+imol
232C     GO TO 541
233C       
234541   continue
235      if(a.le.anu2) go to 19
236      if(kanal.eq.1) go to 16
237      im=im+1
238      anu1=anu2
239      anu2=amin1(nu2,anu2+dnu)
240      go to 20
24116    continue
242      call impanl(tab,anu1,anu2,impr,ymoyd,ymaxd,alf,qq,ibase)
243      anu1=anu2
244      anu2=amin1(nu2,anu2+dnu)
245      do 17 j=1,lk
24617    tab(j)=0
247      do 18 j=1,nmol
248      ymoy(j)=0.
249      ymax(j)=0.
250      ymoyd(j)=0.
251      ymaxd(j)=0.
25218    alf(j)=0.
253      im=im+1
25419    continue
255      if(.not.qq(imol)) go to 100
256      if(kanal.eq.0) go to 20
257      tab(izot)=tab(izot)+1
258      alf(imol)=alf(imol)+v(3)
259      a2=v(2)*(1./cor)
260      ymoyd(imol)=ymoyd(imol)+a2
261Camax1 -> dmax1 ( real 8) BB 19/11/96
262      ymaxd(imol)=dmax1(ymaxd(imol),a2)
263C     ymax(imol)=amax1(ymax(imol),v(2))
264C       
265C     CALCUL DES MIN ET DES MAX DES CLASSES POUR LES INTENSITES ET LES
266C     NIVEAUX DE BASE
267C       
268      if(khist.eq.0.or.khist.eq.1) go to 100
26920    continue
270      if(ival.ne.0.and.kode(iran).ne.ival) go to 100
271      xmax(im+1)=amax1(xmax(im+1),v(4))
272      xmin(im+1)=amin1(xmin(im+1),v(4))
273      go to 100
274200   continue
275      if(kanal.eq.0) go to 24
276      call impanl(tab,anu1,anu2,impr,ymoyd,ymaxd,alf,qq,ibase)
277      if(mode.ne.0.or.modif.ne.oui) go to 900
278c     print *,' avant read rec=1   '
279      read (iuni,rec=1)
280     &nu1,nu2,anu,n203,nbraie,nbmol,iecr,ifin,ll1,ll2,ll3,ll4
281c     print *,' apres read rec=1   '
282c     print *,
283c    &nu1,nu2,anu,n203,nbraie,nbmol,iecr,ifin,ll1,ll2,ll3,ll4
284      ll2=1
285c     print *,' avant write rec=1   '
286c     print *, ifin,ll1
287      write(iuni,rec=1)
288     &nu1,nu2,anu,n203,nbraie,nbmol,iecr,ifin,ll1,ll2,ll3,ll4
289c     print *,' apres write rec=1   '
290c     print *,
291c    &nu1,nu2,anu,n203,nbraie,nbmol,iecr,ifin,ll1,ll2,ll3,ll4     
292c     print *, ifin,ll1
293      ifin=ifin+ll1
294      do 22 j=1,nmol
295      ymoy(j)=ymoyd(j)*cor
29622    ymax(j)=ymaxd(j)*cor
297c     print *, ifin
298c     write(iuni,rec=ifin)
299c    & ntab,(tab(j),j=1,ntab),nmol,(ymoy(j),j=1,nmol),
300c    &(ymax(j),j=1,nmol),(alf(j),j=1,nmol)
301c     print *,
302c    & ntab,(tab(j),j=1,ntab),nmol,(ymoy(j),j=1,nmol),
303c    &(ymax(j),j=1,nmol),(alf(j),j=1,nmol)
304      go to 900
30524    continue
306      impr=0
307      im=0
308      anu1=nu1
309      anu2=nu1+dnu
310      if(dnu.eq.0.) anu2=nu2
311      pas(   2)=(xmax(im+1)-xmin(im+1))/nbclas
312      if(pas(2).eq.0.) pas(2)=1.e+20
313      call pgeisa(nu1,nu2,*900)
314300   continue
315      call lgeisa(v,*400)
316C       
317C     CALCUL DES FREQUENCES DES CLASSES DE HIST
318C       
319      if(a.le.anu2) go to 319
320      call imph  (hist,xmin(im+1),xmax(im+1),pas,anu1,anu2,impr)
321      im=im+1
322      anu1=anu2
323      anu2=amin1(nu2,anu2+dnu)
324      pas(2)=(xmax(im+1)-xmin(im+1))/nbclas
325      if(pas(2).eq.0.) pas(2)=1.e+20
326      do 26 i=1,nbc2
327      hist(i)=0
32826    continue
329319   continue
330      if(ival.ne.0.and.kode(iran).ne.ival) go to 300
331C alog -> dlog (real 8) BB 19/11/96
332      a2=v(2)*(1./cor)
333      kt=incr+dlog10(a2)
334C     kt=incr+alog10(v(2))
335      hist(kt)=hist(kt)+1
336      if(khist.eq.1) go to 300
337      kt=(v(4)-xmin(im+1) )/pas(2) + 1.
338      kt=min0(kt,nbclas)
339      ktk=kt+nbclas
340      hist(ktk)=hist(ktk)+1
341      go to 300
342400   continue
343      if(khist.ne.0) call imph(hist,xmin(im+1),xmax(im+1),pas
344     &,anu1,anu2,impr)
345900   continue
346      return 1
347      end
Note: See TracBrowser for help on using the repository browser.