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    FD218221  25/12/09    21:15:03     12416          
  2.       SUBROUTINE KRIG
  3. C
  4. C     Operateur KRIG (krigeage)
  5. C     Effectue un krigeage ordinaire pour interpoler une fonction
  6. C     sur des points cibles, a partir de valeurs sur des points de mesure
  7.  
  8. C     Typages implicites habituels
  9.       IMPLICIT INTEGER(I-N)
  10.       IMPLICIT REAL*8 (A-H,O-Z)
  11.  
  12. C     Les includes necessaires
  13. -INC PPARAM
  14. -INC CCOPTIO
  15. -INC CCREEL
  16. -INC SMLMOTS
  17. -INC SMLREEL
  18. -INC SMEVOLL
  19. -INC SMELEME
  20. -INC SMCHPOI
  21. -INC SMTABLE
  22.  
  23. C     Un segment pour stocker les donnees
  24. C     NC : nombre de coordonnees des points
  25. C     MM : nombre de points de mesure
  26. C     NN : nombre de points cibles
  27. C     MCOO(i,k) : k-eme coordonnee du i-eme point de mesure
  28. C     MVAL(i,1) : valeur fonction au  i-eme point de mesure
  29. C     CCOO(j,k) : k-eme coordonnee du j-eme point cible
  30. C     DM(i1,i2) : distance entre les points de mesure i1 et i2
  31. C     DC(i,j)   : distance entre le point de mesure i et le point cible j
  32. C     VM(i1,i2) : variance (ou covariance) entre les points de mesure i1 et i2
  33. C                 avec une ligne supplementaire (1 1 1 ... 1 0) pour assurer le krigeage ordinaire
  34. C     VC(i,j)   : Avant appel a DGESV : variance (ou covariance) entre le point mesure i
  35. C                                       et le point cible j
  36. C                                       la derniere ligne est (1 1 1 ... 1) pour le krigeage ordinaire
  37. C                 Apres appel a DGESV : poids d'interpolation pour la i-eme mesure
  38. C                                       pour interpoler au j-eme point cible
  39. C                                       la derniere ligne contient les multiplicateurs de Lagrange
  40. C                                       pour le krigeage ordinaire
  41.       SEGMENT SKRIG
  42.           REAL*8 MCOO(MM,NC)
  43.           REAL*8 MVAL(MM,1)
  44.           REAL*8 CCOO(NN,NC)
  45.           REAL*8 DM(MM,MM)
  46.           REAL*8 DC(MM,NN)
  47.           REAL*8 VM(MM+1,MM+1)
  48.           REAL*8 VC(MM+1,NN)
  49.           REAL*8 WC(MM+1,NN)
  50.           INTEGER IPIV(MM+1)
  51.       ENDSEGMENT
  52.  
  53. C     Quelques SEGMENTs LISTREEL
  54.       POINTEUR MLREE4.MLREEL
  55.  
  56. C     Quelques nombres/objets utiles
  57.       PARAMETER(XUN=1.D0,DEUX=2.D0)
  58.       PARAMETER(XNOR=XPI/180.D0)
  59.       CHARACTER*4 MOT1,COMP
  60.       CHARACTER*11 MDIST
  61.       CHARACTER*8 TYPOBJ
  62.       LOGICAL BVALE
  63.  
  64.  
  65.  
  66. C----------------------------------------------------------------------C
  67. C             ACQUISITION DES DONNEES D'ENTREE UTILISATEUR             C
  68. C----------------------------------------------------------------------C
  69.  
  70. C     La table (ITB)
  71.       CALL LIROBJ('TABLE',ITB,1,IRETOU)
  72.       IF(IERR.NE.0) RETURN
  73.  
  74. C     Donne t'on une distance ?
  75. C     IDIST=1 : Euclidienne (defaut)
  76. C     IDIST=2 : Spherique
  77. C     IDIST=3 : Cylindrique (a faire ...)
  78.       IDIST=1
  79.       TYPOBJ='        '
  80.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'DISTANCE',BIND,IOIND,
  81.      &            TYPOBJ,IVALE,XVALE,MDIST,BVALE,IOVALE)
  82.       IF (IERR.NE.0) GOTO 999
  83.       IF (TYPOBJ.EQ.'MOT     ') THEN
  84.           IF (MDIST.EQ.'EUCLIDIENNE') THEN
  85.               IDIST=1
  86.           ELSEIF(MDIST.EQ.'SPHERIQUE') THEN
  87.               IDIST=2
  88. C             Donne t'on un rayon de sphere ?
  89.               TYPOBJ='        '
  90.               CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'RAYON',BIND,IOIND,
  91.      &                    TYPOBJ,IVALE,XVALE,' ',BVALE,IOVALE)
  92.               IF (TYPOBJ.EQ.'FLOTTANT') THEN
  93.                   RAYON=XVALE
  94.               ELSE
  95.                   RAYON=1.D0
  96.               ENDIF
  97.           ENDIF
  98.       ELSE
  99.           IDIST=1
  100.       ENDIF
  101.  
  102. C     Noms des dimensions (LISTMOTS : MLMOT1)
  103. C     Si distance spherique, on impose la liste 'LONG','LATI'
  104.       IF (IDIST.EQ.2) THEN
  105.           JGN=4
  106.           JGM=2
  107.           SEGINI MLMOT1
  108.           MLMOT1.MOTS(1)='LONG'
  109.           MLMOT1.MOTS(2)='LATI'
  110. C     Pour les autres distances, acquisition obligatoire d'un LISTMOT
  111.       ELSE
  112.           CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'COORDONNEES',BIND,IOIND,
  113.      &                'LISTMOTS',IVALE,XVALE,' ',BVALE,MLMOT1)
  114.           IF (IERR.NE.0) GOTO 999
  115.           SEGACT MLMOT1
  116.       ENDIF
  117.  
  118. C     Mesures (CHPOINT : MCHPO1)
  119.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'MESURES',BIND,IOIND,
  120.      &            'CHPOINT',IVALE,XVALE,' ',BVALE,MCHPO1)
  121.       IF (IERR.NE.0) GOTO 999
  122.       SEGACT MCHPO1
  123. C     Maillage de POI1 support des mesures (IPT1)
  124.       CALL EXTR21(MCHPO1,1,IPT1)
  125. C     Nombre de points de mesure (M)
  126.       M=IPT1.NUM(/2)
  127.  
  128. C     Nom de la composante a kirger (MOT : COMP)
  129.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'COMPOSANTE',BIND,IOIND,
  130.      &            'MOT     ',IVALE,XVALE,COMP,BVALE,IOVALE)
  131.       IF (IERR.NE.0) GOTO 999
  132.  
  133. C     Cibles (CHPOINT : MCHPO2)
  134.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'CIBLES',BIND,IOIND,
  135.      &            'CHPOINT ',IVALE,XVALE,' ',BVALE,MCHPO2)
  136.       IF (IERR.NE.0) GOTO 999
  137.       SEGACT MCHPO2
  138. C     Maillage de POI1 support des cibles (IPT2)
  139.       CALL EXTR21(MCHPO2,1,IPT2)
  140. C     Nombre de points cibles (N)
  141.       N=IPT2.NUM(/2)
  142.  
  143. C     Variogramme ? (IVAR=1) ou covariogramme ? (IVAR=2) (EVOLUTIO : MEVOL1)
  144.       IVAR=0
  145.       TYPOBJ='        '
  146.       CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'VARIOGRAMME',BIND,IOIND,
  147.      &            TYPOBJ,IVALE,XVALE,' ',BVALE,MEVOL1)
  148.       IF (IERR.NE.0) GOTO 999
  149.       IF (TYPOBJ.EQ.'EVOLUTIO') THEN
  150.           IVAR=1
  151.       ELSE
  152.           CALL ACCTAB(ITB,'MOT     ',IIND,XIND,'COVARIOGRAMME',BIND,
  153.      &          IOIND,TYPOBJ,IVALE,XVALE,' ',BVALE,MEVOL1)
  154.           IF (IERR.NE.0) GOTO 999
  155.           IF (TYPOBJ.EQ.'EVOLUTIO') THEN
  156.               IVAR=2
  157.           ELSE
  158.               MOTERR(1:40)='VARIOGRAMME ou COVARIOGRAMME            '
  159.               CALL ERREUR(535)
  160.               RETURN
  161.           ENDIF
  162.       ENDIF
  163. C     Recuperation des 2 LISTREELs du variogramme/covariogramme
  164.       SEGACT MEVOL1
  165.       KEVOL1=MEVOL1.IEVOLL(1)
  166.       SEGACT KEVOL1
  167.       MLREE1=KEVOL1.IPROGX
  168.       MLREE2=KEVOL1.IPROGY
  169.  
  170. C     Rangement des donnees dans un segment SKRIG
  171. C     Nombre de mesures
  172.       MM=M
  173. C     Nombre de cibles
  174.       NN=N
  175. C     Nombre de coordonnees
  176.       NC=MLMOT1.MOTS(/2)
  177.       SEGINI SKRIG
  178. C     Remplissage des coordonnees des points (mesures/cibles)
  179.       DO K=1,NC
  180.           MOT1=MLMOT1.MOTS(K)
  181.           DO I=1,MM
  182. C             Si les mesures sont fournies sous forme de CHPOINT
  183.               IP1=IPT1.NUM(1,I)
  184.               CALL EXTRA9(MCHPO1,IP1,MOT1,IFOUR,.FALSE.,XFLOT,IRET)
  185.               MCOO(I,K)=XFLOT
  186.           ENDDO
  187.           DO I=1,NN
  188. C             Si les cibles sont fournies sous forme de CHPOINT
  189.               IP2=IPT2.NUM(1,I)
  190.               CALL EXTRA9(MCHPO2,IP2,MOT1,IFOUR,.FALSE.,XFLOT,IRET)
  191.               CCOO(I,K)=XFLOT
  192.           ENDDO
  193.       ENDDO
  194. C     Remplissage des valeurs de la fonction a kriger
  195.       DO I=1,MM
  196. C         Si les mesures sont fournies sous forme de CHPOINT
  197.           IP1=IPT1.NUM(1,I)
  198.           CALL EXTRA9(MCHPO1,IP1,COMP,IFOUR,.FALSE.,XFLOT,IRET)
  199.           MVAL(I,1)=XFLOT
  200.       ENDDO
  201.  
  202.  
  203.  
  204. C----------------------------------------------------------------------C
  205. C                         CALCUL DES DISTANCES                         C
  206. C----------------------------------------------------------------------C
  207.  
  208. C     Entre les points de mesure 2 a 2
  209.       DO I=1,MM
  210.           DO J=I,MM
  211. C             Distance euclidienne
  212.               IF (IDIST.EQ.1) THEN
  213.                   DIJ=XZERO
  214.                   DO K=1,NC
  215.                       XI=MCOO(I,K)
  216.                       XJ=MCOO(J,K)
  217.                       DIJ=DIJ+((XI-XJ)**2)
  218.                   ENDDO
  219.                   DIJ=SQRT(DIJ)
  220. C             Distance sur une sphere (formule de haversine pour calculer
  221. C             la distance sur une sphere avec les angles de longitude et latitude)
  222.               ELSEIF (IDIST.EQ.2) THEN
  223.                   XLONGI=XNOR*MCOO(I,1)
  224.                   XLATII=XNOR*MCOO(I,2)
  225.                   XLONGJ=XNOR*MCOO(J,1)
  226.                   XLATIJ=XNOR*MCOO(J,2)
  227.                   DIJ=DEUX*RAYON*ASIN(SQRT(
  228.      &                ((SIN((XLATII-XLATIJ)/DEUX))**2)+
  229.      &                (COS(XLATII)*COS(XLATIJ)*
  230.      &                 ((SIN((XLONGI-XLONGJ)/DEUX))**2)) ))
  231.               ENDIF
  232.               DM(I,J)=DIJ
  233.               DM(J,I)=DIJ
  234.           ENDDO
  235.       ENDDO
  236.  
  237. C     Entre les points cibles et les points de mesure
  238.       DO I=1,MM
  239.           DO J=1,NN
  240. C             Distance euclidienne
  241.               IF (IDIST.EQ.1) THEN
  242.                   DIJ=XZERO
  243.                   DO K=1,NC
  244.                       XI=MCOO(I,K)
  245.                       XJ=CCOO(J,K)
  246.                       DIJ=DIJ+((XI-XJ)**2)
  247.                   ENDDO
  248.                   DIJ=SQRT(DIJ)
  249. C             Distance sur une sphere (formule de haversine pour calculer
  250. C             la distance sur une sphere avec les angles de longitude et latitude)
  251.               ELSEIF (IDIST.EQ.2) THEN
  252.                   XLONGI=XNOR*MCOO(I,1)
  253.                   XLATII=XNOR*MCOO(I,2)
  254.                   XLONGJ=XNOR*CCOO(J,1)
  255.                   XLATIJ=XNOR*CCOO(J,2)
  256.                   DIJ=DEUX*RAYON*ASIN(SQRT(
  257.      &                ((SIN((XLATII-XLATIJ)/DEUX))**2)+
  258.      &                (COS(XLATII)*COS(XLATIJ)*
  259.      &                 ((SIN((XLONGI-XLONGJ)/DEUX))**2))))
  260.               ENDIF
  261.               DC(I,J)=DIJ
  262.           ENDDO
  263.       ENDDO
  264.  
  265.  
  266.  
  267. C----------------------------------------------------------------------C
  268. C                      CALCUL DES (CO)VARIANCES :                      C
  269. C                  MATRICE ET VECTEURS SECOND MEMBRES                  C
  270. C----------------------------------------------------------------------C
  271.  
  272. C     Entre les points de mesure 2 a 2
  273.       DO I=1,MM
  274.           DO J=I,MM
  275.               DIJ=DM(I,J)
  276. C             Interpolation depuis l'EVOLUTIO du variogramme/covariogramme
  277.               CALL INTER5(DIJ,MLREE1,MLREE2,VIJ,0,0,1,IRET)
  278.               IF (IRET.NE.1) GOTO 999
  279.               VM(I,J)=VIJ
  280.               VM(J,I)=VIJ
  281.           ENDDO
  282.       ENDDO
  283. C     Ajout d'une ligne pour imposer le krigeage ordinaire (somme des inconnues = 1)
  284.       DO J=1,MM
  285.           VM(MM+1,J)=XUN
  286.           VM(J,MM+1)=XUN
  287.       ENDDO
  288.       VM(MM+1,MM+1)=XZERO
  289.  
  290. C     Entre les points cibles et les points de mesure
  291.       DO I=1,MM
  292.           DO J=1,NN
  293.               DIJ=DC(I,J)
  294. C             Interpolation depuis l'EVOLUTIO du variogramme/covariogramme
  295.               CALL INTER5(DIJ,MLREE1,MLREE2,VIJ,0,0,1,IRET)
  296.               IF (IRET.NE.1) GOTO 999
  297.               VC(I,J)=VIJ
  298.               WC(I,J)=VIJ
  299.           ENDDO
  300.       ENDDO
  301. C     Ajout d'une ligne pour imposer le krigeage ordinaire (somme des inconnues = 1)
  302.       DO J=1,NN
  303.           VC(MM+1,J)=XUN
  304.           WC(MM+1,J)=XUN
  305.       ENDDO
  306.  
  307.  
  308.  
  309. C----------------------------------------------------------------------C
  310. C                    RESOLUTION DU SYSTEME LINEAIRE                    C
  311. C----------------------------------------------------------------------C
  312.  
  313. C     Resolution du systeme VM * WC = VC a l'aide du solveur de LAPACK
  314.       CALL DGESV(MM+1,NN,VM,MM+1,IPIV,WC,MM+1,INFO)
  315. C     Sortie en cas d'erreur
  316.       IF (INFO.NE.0) THEN
  317.           PRINT*,'ERROR IN DGESV SOLVER, INFO=',INFO
  318.           PRINT*,'CHECK THE LINEAR SYSTEM'
  319.           CALL ERREUR(26)
  320.           RETURN
  321.       ENDIF
  322.  
  323.  
  324.  
  325. C----------------------------------------------------------------------C
  326. C         INTERPOLATION LINEAIRE DES VALEURS AUX POINTS CIBLES         C
  327. C----------------------------------------------------------------------C
  328.  
  329. C     Estimation de la fonction
  330.       JG=NN
  331.       SEGINI MLREE3
  332.       DO J=1,NN
  333.           XVAL=XZERO
  334.           DO I=1,MM
  335.               XVAL=XVAL+(WC(I,J)*MVAL(I,1))
  336.           ENDDO
  337.           MLREE3.PROG(J)=XVAL
  338.       ENDDO
  339.  
  340. C     Variance d'estimation
  341.       JG=NN
  342.       SEGINI MLREE4
  343.       DO J=1,NN
  344. C         Si on a utilise les variances, produit scalaire simple
  345.           IF (IVAR.EQ.1) THEN
  346.               XVAL=XZERO
  347.               DO I=1,MM+1
  348.                     XVAL=XVAL+(WC(I,J)*VC(I,J))
  349.               ENDDO
  350. C         Si on a utilise les covariances, petite soustraction a faire
  351.           ELSEIF (IVAR.EQ.2) THEN
  352.               CALL INTER5(XZERO,MLREE1,MLREE2,XSIG,0,0,1,IRET)
  353.               IF (IRET.NE.1) GOTO 999
  354.               XVAL=XSIG
  355.               DO I=1,MM+1
  356.                     XVAL=XVAL-(WC(I,J)*VC(I,J))
  357.               ENDDO
  358.           ENDIF
  359.           MLREE4.PROG(J)=XVAL
  360.       ENDDO
  361.  
  362.  
  363.  
  364. C----------------------------------------------------------------------C
  365. C                        ECRITURE DES RESULTATS                        C
  366. C----------------------------------------------------------------------C
  367.  
  368. C     Resultats sous forme de CHPOINT
  369. C     Variance d'estimation
  370.       CALL ECROBJ('LISTREEL',MLREE4)
  371.       CALL ECRCHA(COMP)
  372.       CALL ECROBJ('MAILLAGE',IPT2)
  373.       CALL MANUCH
  374.       CALL LIROBJ('CHPOINT',MCHPO4,1,IRETOU)
  375.       SEGACT MCHPO4*MOD
  376.       MCHPO4.JATTRI(1)=MCHPO2.JATTRI(1)
  377.       CALL ECROBJ('CHPOINT',MCHPO4)
  378. C     Estimation de la fonction
  379.       CALL ECROBJ('LISTREEL',MLREE3)
  380.       CALL ECRCHA(COMP)
  381.       CALL ECROBJ('MAILLAGE',IPT2)
  382.       CALL MANUCH
  383.       CALL LIROBJ('CHPOINT',MCHPO3,1,IRETOU)
  384.       SEGACT MCHPO3*MOD
  385.       MCHPO3.JATTRI(1)=MCHPO2.JATTRI(1)
  386.       CALL ECROBJ('CHPOINT',MCHPO3)
  387.  
  388. C     On fait le menage avant de partir
  389.       SEGDES MLMOT1
  390.       SEGDES MEVOL1
  391.       SEGDES KEVOL1
  392.       SEGSUP MLREE3,MLREE4
  393.       SEGDES MCHPO1
  394.       SEGDES MCHPO2,MCHPO3,MCHPO4
  395.  
  396. C     Et c'est fini !
  397.       RETURN
  398.  
  399. C     En cas d'erreur
  400. 999   CONTINUE
  401.       CALL ERREUR(26)
  402.       RETURN
  403.  
  404.       END
  405.  
  406.  

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