C KRIG SOURCE CB215821 25/12/19 21:15:01 12430 SUBROUTINE KRIG C C Gestion de l'operateur KRIG (krigeage) C C Effectue un krigeage ordinaire pour interpoler une fonction C sur des points cibles, a partir de valeurs sur des points de mesure C Les donnees d'entree (mesures / cibles) peuvent etre fournies via C differents objets : C - LISTREEL C - CHPOINT C - MCHAML (a faire ...) C La resolution du systeme lineaire de krigeage aux points cibles est C realisee a l'aide du solveur DGESV (issu de la bibliotheque LAPACK) C Typages implicites habituels IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C Les includes necessaires -INC CCNOYAU -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SMLMOTS -INC SMLREEL -INC SMEVOLL -INC SMELEME -INC SMCHPOI -INC SMTABLE C Un segment pour stocker les donnees C NC : nombre de coordonnees des points C MM : nombre de points de mesure C NN : nombre de points cibles C MCOO(i,k) : k-eme coordonnee du i-eme point de mesure C MVAL(i,1) : valeur fonction au i-eme point de mesure C CCOO(j,k) : k-eme coordonnee du j-eme point cible C DM(i1,i2) : distance entre les points de mesure i1 et i2 C DC(i,j) : distance entre le point de mesure i et le point cible j C VM(i1,i2) : variance (ou covariance) entre les points de mesure i1 et i2 C avec une ligne supplementaire (1 1 1 ... 1 0) pour assurer le krigeage ordinaire C VC(i,j) : Avant appel a DGESV : variance (ou covariance) entre le point mesure i C et le point cible j C la derniere ligne est (1 1 1 ... 1) pour le krigeage ordinaire C Apres appel a DGESV : poids d'interpolation pour la i-eme mesure C pour interpoler au j-eme point cible C la derniere ligne contient les multiplicateurs de Lagrange C pour le krigeage ordinaire SEGMENT SKRIG REAL*8 MCOO(MM,NC) REAL*8 MVAL(MM,1) REAL*8 CCOO(NN,NC) REAL*8 DM(MM,MM) REAL*8 DC(MM,NN) REAL*8 VM(MM+1,MM+1) REAL*8 VC(MM+1,NN) REAL*8 WC(MM+1,NN) INTEGER IPIV(MM+1) ENDSEGMENT C Quelques SEGMENTs LISTREEL POINTEUR MLREE4.MLREEL,MLREE5.MLREEL C Quelques nombres/objets utiles PARAMETER(XUN=1.D0,DEUX=2.D0) PARAMETER(XNOR=XPI/180.D0) CHARACTER*4 MOT4,COMP CHARACTER*11 MDIST CHARACTER*8 MOT8,TYPOBJ LOGICAL BVALE C----------------------------------------------------------------------C C ACQUISITION DES DONNEES D'ENTREE UTILISATEUR C C----------------------------------------------------------------------C C La table (ITB) CALL LIROBJ('TABLE',ITB,1,IRETOU) IF(IERR.NE.0) RETURN C Donne t'on une 'DISTANCE' ? C IDIST=1 : Euclidienne (defaut) C IDIST=2 : Spherique C IDIST=3 : Cylindrique (a faire ...) IDIST=1 TYPOBJ=' ' CALL ACCTAB(ITB,'MOT ',IIND,XIND,'DISTANCE',BIND,IOIND, & TYPOBJ,IVALE,XVALE,MDIST,BVALE,IOVALE) IF (IERR.NE.0) GOTO 999 IF (TYPOBJ.EQ.'MOT ') THEN IF (MDIST.EQ.'EUCLIDIENNE') THEN IDIST=1 ELSEIF(MDIST.EQ.'SPHERIQUE') THEN IDIST=2 C Donne t'on un rayon de sphere ? TYPOBJ=' ' CALL ACCTAB(ITB,'MOT ',IIND,XIND,'RAYON',BIND,IOIND, & TYPOBJ,IVALE,XVALE,' ',BVALE,IOVALE) IF (TYPOBJ.EQ.'FLOTTANT') THEN RAYON=XVALE ELSE RAYON=1.D0 ENDIF ENDIF ELSE IDIST=1 ENDIF C Noms des dimensions ('COORDONNEES') ? C Si distance spherique, on impose la liste 'LONG','LATI' IF (IDIST.EQ.2) THEN JGN=4 JGM=2 SEGINI MLMOT1 MLMOT1.MOTS(1)='LONG' MLMOT1.MOTS(2)='LATI' C Pour les autres distances, acquisition obligatoire d'un LISTMOT ELSE CALL ACCTAB(ITB,'MOT ',IIND,XIND,'COORDONNEES',BIND,IOIND, & 'LISTMOTS',IVALE,XVALE,' ',BVALE,MLMOT1) IF (IERR.NE.0) GOTO 999 SEGACT MLMOT1 ENDIF C Points de 'MESURES' ? C Plusieurs formats possibles : C ITYPM = 1 : TABLE de LISTREELs C ITYPM = 2 : CHPOINT (MCHPO1) ITYPM=0 TYPOBJ=' ' CALL ACCTAB(ITB,'MOT ',IIND,XIND,'MESURES',BIND,IOIND, & TYPOBJ,IVALE,XVALE,' ',BVALE,IOVALE) IF (IERR.NE.0) GOTO 999 C Si mesures dans une TABLE de LISTREELs IF (TYPOBJ.EQ.'TABLE ') THEN ITYPM=1 MTAB1=IOVALE SEGACT MTAB1 TYPOBJ=MTAB1.MTABTI(1) I11=MTAB1.MTABII(1) IDEB=IPCHAR(I11) IFIN=IPCHAR(I11+1) MOT8=ICHARA(IDEB:IFIN-1) C Nombre de points de mesure (M) C pour le trouver, on utilise le premier LISTREEL de la table CALL ACCTAB(MTAB1,'MOT ',IIND,XIND,MOT8,BIND,IOIND, & 'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5) IF (IERR.NE.0) GOTO 999 SEGACT MLREE5 M=MLREE5.PROG(/1) IF (M.LT.1) THEN CALL ERREUR(726) RETURN ENDIF C Si mesures dans un CHPOINT ELSEIF (TYPOBJ.EQ.'CHPOINT ') THEN ITYPM=2 MCHPO1=IOVALE CALL ACTOBJ(TYPOBJ,MCHPO1,1) C Maillage de POI1 support des mesures (IPT1) CALL EXTR21(MCHPO1,1,IPT1) C Nombre de points de mesure (M) M=IPT1.NUM(/2) ELSE MOTERR(1:20)='MESURES LISTREEL ' CALL ERREUR(627) RETURN ENDIF C Nom de la 'COMPOSANTE' a kirger (MOT : COMP) CALL ACCTAB(ITB,'MOT ',IIND,XIND,'COMPOSANTE',BIND,IOIND, & 'MOT ',IVALE,XVALE,COMP,BVALE,IOVALE) IF (IERR.NE.0) GOTO 999 C Points 'CIBLES' ? C Plusieurs formats possibles : C ITYPC = 1 : TABLE de LISTREELs C ITYPC = 2 : CHPOINT (MCHPO2) ITYPC=0 TYPOBJ=' ' CALL ACCTAB(ITB,'MOT ',IIND,XIND,'CIBLES',BIND,IOIND, & TYPOBJ,IVALE,XVALE,' ',BVALE,IOVALE) IF (IERR.NE.0) GOTO 999 C Si cibles dans une TABLE de LISTREELs IF (TYPOBJ.EQ.'TABLE ') THEN ITYPC=1 MTAB2=IOVALE SEGACT MTAB2 TYPOBJ=MTAB2.MTABTI(1) I11=MTAB2.MTABII(1) IDEB=IPCHAR(I11) IFIN=IPCHAR(I11+1) MOT8=ICHARA(IDEB:IFIN-1) C Nombre de points cibles (N) C pour le trouver, on utilise le premier LISTREEL de la table CALL ACCTAB(MTAB2,'MOT ',IIND,XIND,MOT8,BIND,IOIND, & 'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5) IF (IERR.NE.0) GOTO 999 SEGACT MLREE5 N=MLREE5.PROG(/1) IF (N.LT.1) THEN CALL ERREUR(726) RETURN ENDIF C Si cibles dans un CHPOINT ELSEIF (TYPOBJ.EQ.'CHPOINT ') THEN ITYPC=2 MCHPO2=IOVALE CALL ACTOBJ(TYPOBJ,MCHPO2,1) C Maillage de POI1 support des mesures (IPT2) CALL EXTR21(MCHPO2,1,IPT2) C Nombre de points cibles (N) N=IPT2.NUM(/2) ELSE MOTERR(1:20)='CIBLES LISTREEL ' CALL ERREUR(627) RETURN ENDIF C Courbe du 'VARIOGRAMME' ou du 'COVARIOGRAMME' (EVOLUTIO : MEVOL1) ? C IVAR = 1 : variogramme C IVAR = 2 : covariogramme IVAR=0 TYPOBJ=' ' CALL ACCTAB(ITB,'MOT ',IIND,XIND,'VARIOGRAMME',BIND,IOIND, & TYPOBJ,IVALE,XVALE,' ',BVALE,MEVOL1) IF (IERR.NE.0) GOTO 999 IF (TYPOBJ.EQ.'EVOLUTIO') THEN IVAR=1 ELSE CALL ACCTAB(ITB,'MOT ',IIND,XIND,'COVARIOGRAMME',BIND, & IOIND,TYPOBJ,IVALE,XVALE,' ',BVALE,MEVOL1) IF (IERR.NE.0) GOTO 999 IF (TYPOBJ.EQ.'EVOLUTIO') THEN IVAR=2 ELSE MOTERR(1:40)='VARIOGRAMME ou COVARIOGRAMME ' CALL ERREUR(535) RETURN ENDIF ENDIF C Recuperation des 2 LISTREELs du variogramme/covariogramme CALL ACTOBJ('EVOLUTIO',MEVOL1,1) KEVOL1=MEVOL1.IEVOLL(1) MLREE1=KEVOL1.IPROGX MLREE2=KEVOL1.IPROGY C----------------------------------------------------------------------C C RANGEMENT DES DONNEES DANS UN SEGMENT SKRIG C C----------------------------------------------------------------------C C Nombre de mesures MM=M C Nombre de cibles NN=N C Nombre de coordonnees NC=MLMOT1.MOTS(/2) SEGINI SKRIG C Remplissage des coordonnees des points (mesures/cibles) C Boucle sur les noms des coordonnees DO K=1,NC C Coordonnee MOT4 MOT4=MLMOT1.MOTS(K) C Pour les points de mesure C si les donnees des mesures sont dans une TABLE de LISTREELs IF (ITYPM.EQ.1) THEN C Boucle sur les indice de cette table CALL ACCTAB(MTAB1,'MOT ',IIND,XIND,MOT4,BIND,IOIND, & 'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5) IF (IERR.NE.0) GOTO 999 SEGACT MLREE5 C Erreur si les listes n'ont pas la meme taille JG=MLREE5.PROG(/1) IF (JG.NE.MM) THEN CALL ERREUR(217) RETURN ENDIF DO I=1,JG MCOO(I,K)=MLREE5.PROG(I) ENDDO C si les donnees sont dans un CHPOINT ELSEIF (ITYPM.EQ.2) THEN C Boucle sur les points de mesure DO I=1,MM IP1=IPT1.NUM(1,I) CALL EXTRA9(MCHPO1,IP1,MOT4,IFOUR,.FALSE.,XFLOT,IRET) MCOO(I,K)=XFLOT ENDDO ENDIF C Pour les points cibles C si les donnees des cibles sont dans une TABLE de LISTREELs IF (ITYPC.EQ.1) THEN C Boucle sur les indice de cette table CALL ACCTAB(MTAB2,'MOT ',IIND,XIND,MOT4,BIND,IOIND, & 'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5) IF (IERR.NE.0) GOTO 999 SEGACT MLREE5 C Erreur si les listes n'ont pas la meme taille JG=MLREE5.PROG(/1) IF (JG.NE.NN) THEN CALL ERREUR(217) RETURN ENDIF DO I=1,JG CCOO(I,K)=MLREE5.PROG(I) ENDDO C si les donnees sont dans un CHPOINT ELSEIF (ITYPC.EQ.2) THEN C Boucle sur les points de mesure DO I=1,NN IP2=IPT2.NUM(1,I) CALL EXTRA9(MCHPO2,IP2,MOT4,IFOUR,.FALSE.,XFLOT,IRET) CCOO(I,K)=XFLOT ENDDO ENDIF ENDDO C Remplissage des valeurs de la fonction a kriger C Si les donnees des mesures sont dans une TABLE de LISTREELs IF (ITYPM.EQ.1) THEN CALL ACCTAB(MTAB1,'MOT ',IIND,XIND,COMP,BIND,IOIND, & 'LISTREEL',IVALE,XVALE,' ',BVALE,MLREE5) IF (IERR.NE.0) GOTO 999 SEGACT MLREE5 C Erreur si les listes n'ont pas la meme taille JG=MLREE5.PROG(/1) IF (JG.NE.MM) THEN CALL ERREUR(217) RETURN ENDIF DO I=1,JG MVAL(I,1)=MLREE5.PROG(I) ENDDO C Si les donnees sont dans un CHPOINT ELSEIF (ITYPM.EQ.2) THEN C Boucle sur les points de mesure DO I=1,MM IP1=IPT1.NUM(1,I) CALL EXTRA9(MCHPO1,IP1,COMP,IFOUR,.FALSE.,XFLOT,IRET) MVAL(I,1)=XFLOT ENDDO ENDIF C----------------------------------------------------------------------C C CALCUL DES DISTANCES C C----------------------------------------------------------------------C C Entre les points de mesure 2 a 2 DO I=1,MM DO J=I,MM C Distance euclidienne IF (IDIST.EQ.1) THEN DIJ=XZERO DO K=1,NC XI=MCOO(I,K) XJ=MCOO(J,K) DIJ=DIJ+((XI-XJ)**2) ENDDO DIJ=SQRT(DIJ) C Distance sur une sphere (formule de haversine pour calculer C la distance sur une sphere avec les angles de longitude et latitude) ELSEIF (IDIST.EQ.2) THEN XLONGI=XNOR*MCOO(I,1) XLATII=XNOR*MCOO(I,2) XLONGJ=XNOR*MCOO(J,1) XLATIJ=XNOR*MCOO(J,2) DIJ=DEUX*RAYON*ASIN(SQRT( & ((SIN((XLATII-XLATIJ)/DEUX))**2)+ & (COS(XLATII)*COS(XLATIJ)* & ((SIN((XLONGI-XLONGJ)/DEUX))**2)) )) ENDIF DM(I,J)=DIJ DM(J,I)=DIJ ENDDO ENDDO C Entre les points cibles et les points de mesure DO I=1,MM DO J=1,NN C Distance euclidienne IF (IDIST.EQ.1) THEN DIJ=XZERO DO K=1,NC XI=MCOO(I,K) XJ=CCOO(J,K) DIJ=DIJ+((XI-XJ)**2) ENDDO DIJ=SQRT(DIJ) C Distance sur une sphere (formule de haversine pour calculer C la distance sur une sphere avec les angles de longitude et latitude) ELSEIF (IDIST.EQ.2) THEN XLONGI=XNOR*MCOO(I,1) XLATII=XNOR*MCOO(I,2) XLONGJ=XNOR*CCOO(J,1) XLATIJ=XNOR*CCOO(J,2) DIJ=DEUX*RAYON*ASIN(SQRT( & ((SIN((XLATII-XLATIJ)/DEUX))**2)+ & (COS(XLATII)*COS(XLATIJ)* & ((SIN((XLONGI-XLONGJ)/DEUX))**2)))) ENDIF DC(I,J)=DIJ ENDDO ENDDO C----------------------------------------------------------------------C C CALCUL DES (CO)VARIANCES : C C MATRICE ET VECTEURS SECOND MEMBRES C C----------------------------------------------------------------------C C Entre les points de mesure 2 a 2 DO I=1,MM DO J=I,MM DIJ=DM(I,J) C Interpolation depuis l'EVOLUTIO du variogramme/covariogramme CALL INTER5(DIJ,MLREE1,MLREE2,VIJ,0,0,1,IRET) IF (IRET.NE.1) GOTO 999 VM(I,J)=VIJ VM(J,I)=VIJ ENDDO ENDDO C Ajout d'une ligne pour imposer le krigeage ordinaire (somme des inconnues = 1) DO J=1,MM VM(MM+1,J)=XUN VM(J,MM+1)=XUN ENDDO VM(MM+1,MM+1)=XZERO C Entre les points cibles et les points de mesure DO I=1,MM DO J=1,NN DIJ=DC(I,J) C Interpolation depuis l'EVOLUTIO du variogramme/covariogramme CALL INTER5(DIJ,MLREE1,MLREE2,VIJ,0,0,1,IRET) IF (IRET.NE.1) GOTO 999 VC(I,J)=VIJ WC(I,J)=VIJ ENDDO ENDDO C Ajout d'une ligne pour imposer le krigeage ordinaire (somme des inconnues = 1) DO J=1,NN VC(MM+1,J)=XUN WC(MM+1,J)=XUN ENDDO C----------------------------------------------------------------------C C RESOLUTION DU SYSTEME LINEAIRE C C----------------------------------------------------------------------C C Resolution du systeme VM * WC = VC a l'aide du solveur de LAPACK CALL DGESV(MM+1,NN,VM,MM+1,IPIV,WC,MM+1,INFO) C Sortie en cas d'erreur IF (INFO.NE.0) THEN PRINT*,'ERROR IN DGESV SOLVER, INFO=',INFO PRINT*,'CHECK THE LINEAR SYSTEM' CALL ERREUR(26) RETURN ENDIF C----------------------------------------------------------------------C C INTERPOLATION LINEAIRE DES VALEURS AUX POINTS CIBLES C C----------------------------------------------------------------------C C Estimation de la fonction JG=NN SEGINI MLREE3 DO J=1,NN XVAL=XZERO DO I=1,MM XVAL=XVAL+(WC(I,J)*MVAL(I,1)) ENDDO MLREE3.PROG(J)=XVAL ENDDO C Variance d'estimation JG=NN SEGINI MLREE4 DO J=1,NN C Si on a utilise les variances, produit scalaire simple IF (IVAR.EQ.1) THEN XVAL=XZERO DO I=1,MM+1 XVAL=XVAL+(WC(I,J)*VC(I,J)) ENDDO C Si on a utilise les covariances, petite soustraction a faire ELSEIF (IVAR.EQ.2) THEN CALL INTER5(XZERO,MLREE1,MLREE2,XSIG,0,0,1,IRET) IF (IRET.NE.1) GOTO 999 XVAL=XSIG DO I=1,MM+1 XVAL=XVAL-(WC(I,J)*VC(I,J)) ENDDO ENDIF MLREE4.PROG(J)=XVAL ENDDO C----------------------------------------------------------------------C C ECRITURE DES RESULTATS C C----------------------------------------------------------------------C C Si les donnees des cibles sont dans une TABLE de LISTREELs IF (ITYPC.EQ.1) THEN C Variance d'estimation CALL ACTOBJ('LISTREEL',MLREE4,1) CALL ECROBJ('LISTREEL',MLREE4) C Estimation de la fonction CALL ACTOBJ('LISTREEL',MLREE3,1) CALL ECROBJ('LISTREEL',MLREE3) C Si les donnees sont dans un CHPOINT ELSEIF (ITYPC.EQ.2) THEN C Variance d'estimation CALL ECROBJ('LISTREEL',MLREE4) CALL ECRCHA(COMP) CALL ECROBJ('MAILLAGE',IPT2) CALL MANUCH CALL LIROBJ('CHPOINT',MCHPO4,1,IRETOU) SEGACT MCHPO4*MOD MCHPO4.JATTRI(1)=MCHPO2.JATTRI(1) CALL ECROBJ('CHPOINT',MCHPO4) C Estimation de la fonction CALL ECROBJ('LISTREEL',MLREE3) CALL ECRCHA(COMP) CALL ECROBJ('MAILLAGE',IPT2) CALL MANUCH CALL LIROBJ('CHPOINT',MCHPO3,1,IRETOU) SEGACT MCHPO3*MOD MCHPO3.JATTRI(1)=MCHPO2.JATTRI(1) CALL ACTOBJ('CHPOINT',MCHPO3,1) CALL ECROBJ('CHPOINT',MCHPO3) ENDIF C On fait le menage avant de partir IF (ITYPC.EQ.2) THEN SEGSUP MLREE3,MLREE4 ENDIF C Et c'est fini ! RETURN C En cas d'erreur 999 CONTINUE CALL ERREUR(26) RETURN END