Télécharger la plaquette Cast3M


fr en

Télécharger krig.eso

Retour à la liste

Numérotation des lignes : 

krig

  1. C KRIG      SOURCE    CB215821  25/12/19    21:15:01     12430          
  2.       SUBROUTINE KRIG
  3. C
  4. C     Gestion de l'operateur KRIG (krigeage)
  5. C
  6. C     Effectue un krigeage ordinaire pour interpoler une fonction
  7. C     sur des points cibles, a partir de valeurs sur des points de mesure
  8.  
  9. C     Les donnees d'entree (mesures / cibles) peuvent etre fournies via
  10. C     differents objets :
  11. C     - LISTREEL
  12. C     - CHPOINT
  13. C     - MCHAML (a faire ...)
  14.  
  15. C     La resolution du systeme lineaire de krigeage aux points cibles est
  16. C     realisee a l'aide du solveur DGESV (issu de la bibliotheque LAPACK)
  17.  
  18.  
  19.  
  20. C     Typages implicites habituels
  21.       IMPLICIT INTEGER(I-N)
  22.       IMPLICIT REAL*8 (A-H,O-Z)
  23.  
  24. C     Les includes necessaires
  25. -INC CCNOYAU
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC CCREEL
  29. -INC SMLMOTS
  30. -INC SMLREEL
  31. -INC SMEVOLL
  32. -INC SMELEME
  33. -INC SMCHPOI
  34. -INC SMTABLE
  35.  
  36. C     Un segment pour stocker les donnees
  37. C     NC : nombre de coordonnees des points
  38. C     MM : nombre de points de mesure
  39. C     NN : nombre de points cibles
  40. C     MCOO(i,k) : k-eme coordonnee du i-eme point de mesure
  41. C     MVAL(i,1) : valeur fonction au  i-eme point de mesure
  42. C     CCOO(j,k) : k-eme coordonnee du j-eme point cible
  43. C     DM(i1,i2) : distance entre les points de mesure i1 et i2
  44. C     DC(i,j)   : distance entre le point de mesure i et le point cible j
  45. C     VM(i1,i2) : variance (ou covariance) entre les points de mesure i1 et i2
  46. C                 avec une ligne supplementaire (1 1 1 ... 1 0) pour assurer le krigeage ordinaire
  47. C     VC(i,j)   : Avant appel a DGESV : variance (ou covariance) entre le point mesure i
  48. C                                       et le point cible j
  49. C                                       la derniere ligne est (1 1 1 ... 1) pour le krigeage ordinaire
  50. C                 Apres appel a DGESV : poids d'interpolation pour la i-eme mesure
  51. C                                       pour interpoler au j-eme point cible
  52. C                                       la derniere ligne contient les multiplicateurs de Lagrange
  53. C                                       pour le krigeage ordinaire
  54.       SEGMENT SKRIG
  55.           REAL*8 MCOO(MM,NC)
  56.           REAL*8 MVAL(MM,1)
  57.           REAL*8 CCOO(NN,NC)
  58.           REAL*8 DM(MM,MM)
  59.           REAL*8 DC(MM,NN)
  60.           REAL*8 VM(MM+1,MM+1)
  61.           REAL*8 VC(MM+1,NN)
  62.           REAL*8 WC(MM+1,NN)
  63.           INTEGER IPIV(MM+1)
  64.       ENDSEGMENT
  65.  
  66. C     Quelques SEGMENTs LISTREEL
  67.       POINTEUR MLREE4.MLREEL,MLREE5.MLREEL
  68.  
  69. C     Quelques nombres/objets utiles
  70.       PARAMETER(XUN=1.D0,DEUX=2.D0)
  71.       PARAMETER(XNOR=XPI/180.D0)
  72.       CHARACTER*4 MOT4,COMP
  73.       CHARACTER*11 MDIST
  74.       CHARACTER*8 MOT8,TYPOBJ
  75.       LOGICAL BVALE
  76.  
  77.  
  78.  
  79. C----------------------------------------------------------------------C
  80. C             ACQUISITION DES DONNEES D'ENTREE UTILISATEUR             C
  81. C----------------------------------------------------------------------C
  82.  
  83. C     La table (ITB)
  84.       CALL LIROBJ('TABLE',ITB,1,IRETOU)
  85.       IF(IERR.NE.0) RETURN
  86.  
  87. C     Donne t'on une 'DISTANCE' ?
  88. C     IDIST=1 : Euclidienne (defaut)
  89. C     IDIST=2 : Spherique
  90. C     IDIST=3 : Cylindrique (a faire ...)
  91.       IDIST=1
  92.       TYPOBJ='        '
  93.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'DISTANCE',BIND,IOIND,
  94.      &            TYPOBJ,IVALE,XVALE,MDIST,BVALE,IOVALE)
  95.       IF (IERR.NE.0) GOTO 999
  96.       IF (TYPOBJ.EQ.'MOT     ') THEN
  97.           IF (MDIST.EQ.'EUCLIDIENNE') THEN
  98.               IDIST=1
  99.           ELSEIF(MDIST.EQ.'SPHERIQUE') THEN
  100.               IDIST=2
  101. C             Donne t'on un rayon de sphere ?
  102.               TYPOBJ='        '
  103.               CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'RAYON',BIND,IOIND,
  104.      &                    TYPOBJ,IVALE,XVALE,' ',BVALE,IOVALE)
  105.               IF (TYPOBJ.EQ.'FLOTTANT') THEN
  106.                   RAYON=XVALE
  107.               ELSE
  108.                   RAYON=1.D0
  109.               ENDIF
  110.           ENDIF
  111.       ELSE
  112.           IDIST=1
  113.       ENDIF
  114.  
  115. C     Noms des dimensions ('COORDONNEES') ?
  116. C     Si distance spherique, on impose la liste 'LONG','LATI'
  117.       IF (IDIST.EQ.2) THEN
  118.           JGN=4
  119.           JGM=2
  120.           SEGINI MLMOT1
  121.           MLMOT1.MOTS(1)='LONG'
  122.           MLMOT1.MOTS(2)='LATI'
  123. C     Pour les autres distances, acquisition obligatoire d'un LISTMOT
  124.       ELSE
  125.           CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'COORDONNEES',BIND,IOIND,
  126.      &                'LISTMOTS',IVALE,XVALE,' ',BVALE,MLMOT1)
  127.           IF (IERR.NE.0) GOTO 999
  128.           SEGACT MLMOT1
  129.       ENDIF
  130.  
  131. C     Points de 'MESURES' ?
  132. C     Plusieurs formats possibles :
  133. C     ITYPM = 1 : TABLE de LISTREELs
  134. C     ITYPM = 2 : CHPOINT  (MCHPO1)
  135.       ITYPM=0
  136.       TYPOBJ='        '
  137.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'MESURES',BIND,IOIND,
  138.      &            TYPOBJ,IVALE,XVALE,' ',BVALE,IOVALE)
  139.       IF (IERR.NE.0) GOTO 999
  140. C     Si mesures dans une TABLE de LISTREELs
  141.       IF (TYPOBJ.EQ.'TABLE   ') THEN
  142.           ITYPM=1
  143.           MTAB1=IOVALE
  144.           SEGACT MTAB1
  145.           TYPOBJ=MTAB1.MTABTI(1)
  146.           I11=MTAB1.MTABII(1)
  147.           IDEB=IPCHAR(I11)
  148.           IFIN=IPCHAR(I11+1)
  149.           MOT8=ICHARA(IDEB:IFIN-1)
  150. C         Nombre de points de mesure (M)
  151. C         pour le trouver, on utilise le premier LISTREEL de la table
  152.           CALL ACCTAB(MTAB1,'MOT     ',IIND,XIND,MOT8,BIND,IOIND,
  153.      &                'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5)
  154.           IF (IERR.NE.0) GOTO 999
  155.           SEGACT MLREE5
  156.           M=MLREE5.PROG(/1)
  157.           IF (M.LT.1) THEN
  158.               CALL ERREUR(726)
  159.               RETURN
  160.           ENDIF
  161. C     Si mesures dans un CHPOINT
  162.       ELSEIF (TYPOBJ.EQ.'CHPOINT ') THEN
  163.           ITYPM=2
  164.           MCHPO1=IOVALE
  165.           CALL ACTOBJ(TYPOBJ,MCHPO1,1)
  166. C         Maillage de POI1 support des mesures (IPT1)
  167.           CALL EXTR21(MCHPO1,1,IPT1)
  168. C         Nombre de points de mesure (M)
  169.           M=IPT1.NUM(/2)
  170.       ELSE
  171.           MOTERR(1:20)='MESURES    LISTREEL '
  172.           CALL ERREUR(627)
  173.           RETURN
  174.       ENDIF
  175.  
  176. C     Nom de la 'COMPOSANTE' a kirger (MOT : COMP)
  177.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'COMPOSANTE',BIND,IOIND,
  178.      &            'MOT     ',IVALE,XVALE,COMP,BVALE,IOVALE)
  179.       IF (IERR.NE.0) GOTO 999
  180.  
  181. C     Points 'CIBLES' ?
  182. C     Plusieurs formats possibles :
  183. C     ITYPC = 1 : TABLE de LISTREELs
  184. C     ITYPC = 2 : CHPOINT  (MCHPO2)
  185.       ITYPC=0
  186.       TYPOBJ='        '
  187.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'CIBLES',BIND,IOIND,
  188.      &            TYPOBJ,IVALE,XVALE,' ',BVALE,IOVALE)
  189.       IF (IERR.NE.0) GOTO 999
  190. C     Si cibles dans une TABLE de LISTREELs
  191.       IF (TYPOBJ.EQ.'TABLE   ') THEN
  192.           ITYPC=1
  193.           MTAB2=IOVALE
  194.           SEGACT MTAB2
  195.           TYPOBJ=MTAB2.MTABTI(1)
  196.           I11=MTAB2.MTABII(1)
  197.           IDEB=IPCHAR(I11)
  198.           IFIN=IPCHAR(I11+1)
  199.           MOT8=ICHARA(IDEB:IFIN-1)
  200. C         Nombre de points cibles (N)
  201. C         pour le trouver, on utilise le premier LISTREEL de la table
  202.           CALL ACCTAB(MTAB2,'MOT     ',IIND,XIND,MOT8,BIND,IOIND,
  203.      &                'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5)
  204.           IF (IERR.NE.0) GOTO 999
  205.           SEGACT MLREE5
  206.           N=MLREE5.PROG(/1)
  207.           IF (N.LT.1) THEN
  208.               CALL ERREUR(726)
  209.               RETURN
  210.           ENDIF
  211. C     Si cibles dans un CHPOINT
  212.       ELSEIF (TYPOBJ.EQ.'CHPOINT ') THEN
  213.           ITYPC=2
  214.           MCHPO2=IOVALE
  215.           CALL ACTOBJ(TYPOBJ,MCHPO2,1)
  216. C         Maillage de POI1 support des mesures (IPT2)
  217.           CALL EXTR21(MCHPO2,1,IPT2)
  218. C         Nombre de points cibles (N)
  219.           N=IPT2.NUM(/2)
  220.       ELSE
  221.           MOTERR(1:20)='CIBLES     LISTREEL '
  222.           CALL ERREUR(627)
  223.           RETURN
  224.       ENDIF
  225.  
  226. C     Courbe du 'VARIOGRAMME' ou du 'COVARIOGRAMME' (EVOLUTIO : MEVOL1) ?
  227. C     IVAR = 1 : variogramme
  228. C     IVAR = 2 : covariogramme
  229.       IVAR=0
  230.       TYPOBJ='        '
  231.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'VARIOGRAMME',BIND,IOIND,
  232.      &            TYPOBJ,IVALE,XVALE,' ',BVALE,MEVOL1)
  233.       IF (IERR.NE.0) GOTO 999
  234.       IF (TYPOBJ.EQ.'EVOLUTIO') THEN
  235.           IVAR=1
  236.       ELSE
  237.           CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'COVARIOGRAMME',BIND,
  238.      &          IOIND,TYPOBJ,IVALE,XVALE,' ',BVALE,MEVOL1)
  239.           IF (IERR.NE.0) GOTO 999
  240.           IF (TYPOBJ.EQ.'EVOLUTIO') THEN
  241.               IVAR=2
  242.           ELSE
  243.               MOTERR(1:40)='VARIOGRAMME ou COVARIOGRAMME            '
  244.               CALL ERREUR(535)
  245.               RETURN
  246.           ENDIF
  247.       ENDIF
  248. C     Recuperation des 2 LISTREELs du variogramme/covariogramme
  249.       CALL ACTOBJ('EVOLUTIO',MEVOL1,1)
  250.       KEVOL1=MEVOL1.IEVOLL(1)
  251.       MLREE1=KEVOL1.IPROGX
  252.       MLREE2=KEVOL1.IPROGY
  253.  
  254.  
  255.  
  256.  
  257. C----------------------------------------------------------------------C
  258. C             RANGEMENT DES DONNEES DANS UN SEGMENT SKRIG              C
  259. C----------------------------------------------------------------------C
  260.  
  261. C     Nombre de mesures
  262.       MM=M
  263. C     Nombre de cibles
  264.       NN=N
  265. C     Nombre de coordonnees
  266.       NC=MLMOT1.MOTS(/2)
  267.       SEGINI SKRIG
  268.  
  269. C     Remplissage des coordonnees des points (mesures/cibles)
  270. C     Boucle sur les noms des coordonnees
  271.       DO K=1,NC
  272. C         Coordonnee MOT4
  273.           MOT4=MLMOT1.MOTS(K)
  274. C         Pour les points de mesure
  275. C         si les donnees des mesures sont dans une TABLE de LISTREELs
  276.           IF (ITYPM.EQ.1) THEN
  277. C             Boucle sur les indice de cette table
  278.               CALL ACCTAB(MTAB1,'MOT     ',IIND,XIND,MOT4,BIND,IOIND,
  279.      &                    'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5)
  280.               IF (IERR.NE.0) GOTO 999
  281.               SEGACT MLREE5
  282. C             Erreur si les listes n'ont pas la meme taille
  283.               JG=MLREE5.PROG(/1)
  284.               IF (JG.NE.MM) THEN
  285.                   CALL ERREUR(217)
  286.                   RETURN
  287.               ENDIF
  288.               DO I=1,JG
  289.                   MCOO(I,K)=MLREE5.PROG(I)
  290.               ENDDO
  291. C         si les donnees sont dans un CHPOINT
  292.           ELSEIF (ITYPM.EQ.2) THEN
  293. C             Boucle sur les points de mesure
  294.               DO I=1,MM
  295.                   IP1=IPT1.NUM(1,I)
  296.                   CALL EXTRA9(MCHPO1,IP1,MOT4,IFOUR,.FALSE.,XFLOT,IRET)
  297.                   MCOO(I,K)=XFLOT
  298.               ENDDO
  299.           ENDIF
  300. C         Pour les points cibles
  301. C         si les donnees des cibles sont dans une TABLE de LISTREELs
  302.           IF (ITYPC.EQ.1) THEN
  303. C             Boucle sur les indice de cette table
  304.               CALL ACCTAB(MTAB2,'MOT     ',IIND,XIND,MOT4,BIND,IOIND,
  305.      &                    'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5)
  306.               IF (IERR.NE.0) GOTO 999
  307.               SEGACT MLREE5
  308. C             Erreur si les listes n'ont pas la meme taille
  309.               JG=MLREE5.PROG(/1)
  310.               IF (JG.NE.NN) THEN
  311.                   CALL ERREUR(217)
  312.                   RETURN
  313.               ENDIF
  314.               DO I=1,JG
  315.                   CCOO(I,K)=MLREE5.PROG(I)
  316.               ENDDO
  317. C         si les donnees sont dans un CHPOINT
  318.           ELSEIF (ITYPC.EQ.2) THEN
  319. C             Boucle sur les points de mesure
  320.               DO I=1,NN
  321.                   IP2=IPT2.NUM(1,I)
  322.                   CALL EXTRA9(MCHPO2,IP2,MOT4,IFOUR,.FALSE.,XFLOT,IRET)
  323.                   CCOO(I,K)=XFLOT
  324.               ENDDO
  325.           ENDIF
  326.       ENDDO
  327.  
  328. C     Remplissage des valeurs de la fonction a kriger
  329. C     Si les donnees des mesures sont dans une TABLE de LISTREELs
  330.       IF (ITYPM.EQ.1) THEN
  331.           CALL ACCTAB(MTAB1,'MOT     ',IIND,XIND,COMP,BIND,IOIND,
  332.      &                'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5)
  333.           IF (IERR.NE.0) GOTO 999
  334.           SEGACT MLREE5
  335. C         Erreur si les listes n'ont pas la meme taille
  336.           JG=MLREE5.PROG(/1)
  337.           IF (JG.NE.MM) THEN
  338.               CALL ERREUR(217)
  339.               RETURN
  340.           ENDIF
  341.           DO I=1,JG
  342.               MVAL(I,1)=MLREE5.PROG(I)
  343.           ENDDO
  344. C     Si les donnees sont dans un CHPOINT
  345.       ELSEIF (ITYPM.EQ.2) THEN
  346. C         Boucle sur les points de mesure
  347.           DO I=1,MM
  348.               IP1=IPT1.NUM(1,I)
  349.               CALL EXTRA9(MCHPO1,IP1,COMP,IFOUR,.FALSE.,XFLOT,IRET)
  350.               MVAL(I,1)=XFLOT
  351.           ENDDO
  352.       ENDIF
  353.  
  354.  
  355.  
  356. C----------------------------------------------------------------------C
  357. C                         CALCUL DES DISTANCES                         C
  358. C----------------------------------------------------------------------C
  359.  
  360. C     Entre les points de mesure 2 a 2
  361.       DO I=1,MM
  362.           DO J=I,MM
  363. C             Distance euclidienne
  364.               IF (IDIST.EQ.1) THEN
  365.                   DIJ=XZERO
  366.                   DO K=1,NC
  367.                       XI=MCOO(I,K)
  368.                       XJ=MCOO(J,K)
  369.                       DIJ=DIJ+((XI-XJ)**2)
  370.                   ENDDO
  371.                   DIJ=SQRT(DIJ)
  372. C             Distance sur une sphere (formule de haversine pour calculer
  373. C             la distance sur une sphere avec les angles de longitude et latitude)
  374.               ELSEIF (IDIST.EQ.2) THEN
  375.                   XLONGI=XNOR*MCOO(I,1)
  376.                   XLATII=XNOR*MCOO(I,2)
  377.                   XLONGJ=XNOR*MCOO(J,1)
  378.                   XLATIJ=XNOR*MCOO(J,2)
  379.                   DIJ=DEUX*RAYON*ASIN(SQRT(
  380.      &                ((SIN((XLATII-XLATIJ)/DEUX))**2)+
  381.      &                (COS(XLATII)*COS(XLATIJ)*
  382.      &                 ((SIN((XLONGI-XLONGJ)/DEUX))**2)) ))
  383.               ENDIF
  384.               DM(I,J)=DIJ
  385.               DM(J,I)=DIJ
  386.           ENDDO
  387.       ENDDO
  388.  
  389. C     Entre les points cibles et les points de mesure
  390.       DO I=1,MM
  391.           DO J=1,NN
  392. C             Distance euclidienne
  393.               IF (IDIST.EQ.1) THEN
  394.                   DIJ=XZERO
  395.                   DO K=1,NC
  396.                       XI=MCOO(I,K)
  397.                       XJ=CCOO(J,K)
  398.                       DIJ=DIJ+((XI-XJ)**2)
  399.                   ENDDO
  400.                   DIJ=SQRT(DIJ)
  401. C             Distance sur une sphere (formule de haversine pour calculer
  402. C             la distance sur une sphere avec les angles de longitude et latitude)
  403.               ELSEIF (IDIST.EQ.2) THEN
  404.                   XLONGI=XNOR*MCOO(I,1)
  405.                   XLATII=XNOR*MCOO(I,2)
  406.                   XLONGJ=XNOR*CCOO(J,1)
  407.                   XLATIJ=XNOR*CCOO(J,2)
  408.                   DIJ=DEUX*RAYON*ASIN(SQRT(
  409.      &                ((SIN((XLATII-XLATIJ)/DEUX))**2)+
  410.      &                (COS(XLATII)*COS(XLATIJ)*
  411.      &                 ((SIN((XLONGI-XLONGJ)/DEUX))**2))))
  412.               ENDIF
  413.               DC(I,J)=DIJ
  414.           ENDDO
  415.       ENDDO
  416.  
  417.  
  418.  
  419. C----------------------------------------------------------------------C
  420. C                      CALCUL DES (CO)VARIANCES :                      C
  421. C                  MATRICE ET VECTEURS SECOND MEMBRES                  C
  422. C----------------------------------------------------------------------C
  423.  
  424. C     Entre les points de mesure 2 a 2
  425.       DO I=1,MM
  426.           DO J=I,MM
  427.               DIJ=DM(I,J)
  428. C             Interpolation depuis l'EVOLUTIO du variogramme/covariogramme
  429.               CALL INTER5(DIJ,MLREE1,MLREE2,VIJ,0,0,1,IRET)
  430.               IF (IRET.NE.1) GOTO 999
  431.               VM(I,J)=VIJ
  432.               VM(J,I)=VIJ
  433.           ENDDO
  434.       ENDDO
  435. C     Ajout d'une ligne pour imposer le krigeage ordinaire (somme des inconnues = 1)
  436.       DO J=1,MM
  437.           VM(MM+1,J)=XUN
  438.           VM(J,MM+1)=XUN
  439.       ENDDO
  440.       VM(MM+1,MM+1)=XZERO
  441.  
  442. C     Entre les points cibles et les points de mesure
  443.       DO I=1,MM
  444.           DO J=1,NN
  445.               DIJ=DC(I,J)
  446. C             Interpolation depuis l'EVOLUTIO du variogramme/covariogramme
  447.               CALL INTER5(DIJ,MLREE1,MLREE2,VIJ,0,0,1,IRET)
  448.               IF (IRET.NE.1) GOTO 999
  449.               VC(I,J)=VIJ
  450.               WC(I,J)=VIJ
  451.           ENDDO
  452.       ENDDO
  453. C     Ajout d'une ligne pour imposer le krigeage ordinaire (somme des inconnues = 1)
  454.       DO J=1,NN
  455.           VC(MM+1,J)=XUN
  456.           WC(MM+1,J)=XUN
  457.       ENDDO
  458.  
  459.  
  460.  
  461. C----------------------------------------------------------------------C
  462. C                    RESOLUTION DU SYSTEME LINEAIRE                    C
  463. C----------------------------------------------------------------------C
  464.  
  465. C     Resolution du systeme VM * WC = VC a l'aide du solveur de LAPACK
  466.       CALL DGESV(MM+1,NN,VM,MM+1,IPIV,WC,MM+1,INFO)
  467. C     Sortie en cas d'erreur
  468.       IF (INFO.NE.0) THEN
  469.           PRINT*,'ERROR IN DGESV SOLVER, INFO=',INFO
  470.           PRINT*,'CHECK THE LINEAR SYSTEM'
  471.           CALL ERREUR(26)
  472.           RETURN
  473.       ENDIF
  474.  
  475.  
  476.  
  477. C----------------------------------------------------------------------C
  478. C         INTERPOLATION LINEAIRE DES VALEURS AUX POINTS CIBLES         C
  479. C----------------------------------------------------------------------C
  480.  
  481. C     Estimation de la fonction
  482.       JG=NN
  483.       SEGINI MLREE3
  484.       DO J=1,NN
  485.           XVAL=XZERO
  486.           DO I=1,MM
  487.               XVAL=XVAL+(WC(I,J)*MVAL(I,1))
  488.           ENDDO
  489.           MLREE3.PROG(J)=XVAL
  490.       ENDDO
  491.  
  492. C     Variance d'estimation
  493.       JG=NN
  494.       SEGINI MLREE4
  495.       DO J=1,NN
  496. C         Si on a utilise les variances, produit scalaire simple
  497.           IF (IVAR.EQ.1) THEN
  498.               XVAL=XZERO
  499.               DO I=1,MM+1
  500.                     XVAL=XVAL+(WC(I,J)*VC(I,J))
  501.               ENDDO
  502. C         Si on a utilise les covariances, petite soustraction a faire
  503.           ELSEIF (IVAR.EQ.2) THEN
  504.               CALL INTER5(XZERO,MLREE1,MLREE2,XSIG,0,0,1,IRET)
  505.               IF (IRET.NE.1) GOTO 999
  506.               XVAL=XSIG
  507.               DO I=1,MM+1
  508.                     XVAL=XVAL-(WC(I,J)*VC(I,J))
  509.               ENDDO
  510.           ENDIF
  511.           MLREE4.PROG(J)=XVAL
  512.       ENDDO
  513.  
  514.  
  515.  
  516. C----------------------------------------------------------------------C
  517. C                        ECRITURE DES RESULTATS                        C
  518. C----------------------------------------------------------------------C
  519.  
  520. C     Si les donnees des cibles sont dans une TABLE de LISTREELs
  521.       IF (ITYPC.EQ.1) THEN
  522. C         Variance d'estimation
  523.           CALL ACTOBJ('LISTREEL',MLREE4,1)
  524.           CALL ECROBJ('LISTREEL',MLREE4)
  525. C         Estimation de la fonction
  526.           CALL ACTOBJ('LISTREEL',MLREE3,1)
  527.           CALL ECROBJ('LISTREEL',MLREE3)
  528. C     Si les donnees sont dans un CHPOINT
  529.       ELSEIF (ITYPC.EQ.2) THEN
  530. C         Variance d'estimation
  531.           CALL ECROBJ('LISTREEL',MLREE4)
  532.           CALL ECRCHA(COMP)
  533.           CALL ECROBJ('MAILLAGE',IPT2)
  534.           CALL MANUCH
  535.           CALL LIROBJ('CHPOINT',MCHPO4,1,IRETOU)
  536.           SEGACT MCHPO4*MOD
  537.           MCHPO4.JATTRI(1)=MCHPO2.JATTRI(1)
  538.           CALL ECROBJ('CHPOINT',MCHPO4)
  539. C         Estimation de la fonction
  540.           CALL ECROBJ('LISTREEL',MLREE3)
  541.           CALL ECRCHA(COMP)
  542.           CALL ECROBJ('MAILLAGE',IPT2)
  543.           CALL MANUCH
  544.           CALL LIROBJ('CHPOINT',MCHPO3,1,IRETOU)
  545.           SEGACT MCHPO3*MOD
  546.           MCHPO3.JATTRI(1)=MCHPO2.JATTRI(1)
  547.           CALL ACTOBJ('CHPOINT',MCHPO3,1)
  548.           CALL ECROBJ('CHPOINT',MCHPO3)
  549.       ENDIF
  550.  
  551. C     On fait le menage avant de partir
  552.       IF (ITYPC.EQ.2) THEN
  553.           SEGSUP MLREE3,MLREE4
  554.       ENDIF
  555.  
  556. C     Et c'est fini !
  557.       RETURN
  558.  
  559. C     En cas d'erreur
  560. 999   CONTINUE
  561.       CALL ERREUR(26)
  562.       RETURN
  563.  
  564.       END
  565.  
  566.  
  567.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales