|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######spl 00002 MODULE MODE_SPLINES 00003 00004 USE MODI_ABOR1_SFX 00005 USE MODD_SPLINES 00006 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00007 USE PARKIND1 ,ONLY : JPRB 00008 ! 00009 INTERFACE SSCOPY 00010 MODULE PROCEDURE SSCOPY_1 00011 MODULE PROCEDURE SSCOPY_2 00012 END INTERFACE 00013 INTERFACE MXIDML 00014 MODULE PROCEDURE MXIDML 00015 END INTERFACE 00016 INTERFACE MTXAXM 00017 MODULE PROCEDURE MTXAXM 00018 END INTERFACE 00019 INTERFACE MXMSPL 00020 MODULE PROCEDURE MXMSPL_1 00021 MODULE PROCEDURE MXMSPL_2 00022 MODULE PROCEDURE MXMSPL_3 00023 END INTERFACE 00024 INTERFACE SMXINV 00025 MODULE PROCEDURE SMXINV 00026 END INTERFACE 00027 !génération de matrices 00028 INTERFACE SPLIE 00029 MODULE PROCEDURE SPLIE 00030 END INTERFACE 00031 INTERFACE SPLK 00032 MODULE PROCEDURE SPLK_1 00033 MODULE PROCEDURE SPLK_2 00034 END INTERFACE 00035 INTERFACE SPLBFIN 00036 MODULE PROCEDURE SPLBFIN 00037 END INTERFACE 00038 INTERFACE SPLT 00039 MODULE PROCEDURE SPLT_1 00040 MODULE PROCEDURE SPLT_2 00041 END INTERFACE 00042 INTERFACE SPLE 00043 MODULE PROCEDURE SPLE 00044 END INTERFACE 00045 INTERFACE SPLB2E1 00046 MODULE PROCEDURE SPLB2E1 00047 END INTERFACE 00048 INTERFACE SPLB2E 00049 MODULE PROCEDURE SPLB2E 00050 END INTERFACE 00051 INTERFACE TRED2 00052 MODULE PROCEDURE TRED2 00053 END INTERFACE 00054 INTERFACE TQL2_2 00055 MODULE PROCEDURE TQL2_2 00056 END INTERFACE 00057 INTERFACE EISRS1 00058 MODULE PROCEDURE EISRS1 00059 END INTERFACE 00060 INTERFACE SPLV 00061 MODULE PROCEDURE SPLV 00062 END INTERFACE 00063 INTERFACE SPLS2VI 00064 MODULE PROCEDURE SPLS2VI 00065 END INTERFACE 00066 INTERFACE SPLBVM 00067 MODULE PROCEDURE SPLBVM 00068 END INTERFACE 00069 INTERFACE SPLS2V 00070 MODULE PROCEDURE SPLS2V 00071 END INTERFACE 00072 INTERFACE SPLDS2V 00073 MODULE PROCEDURE SPLDS2V 00074 END INTERFACE 00075 INTERFACE SPLPS2V 00076 MODULE PROCEDURE SPLPS2V 00077 END INTERFACE 00078 INTERFACE SPLRI 00079 MODULE PROCEDURE SPLRI 00080 END INTERFACE 00081 INTERFACE SPLRS 00082 MODULE PROCEDURE SPLRS 00083 END INTERFACE 00084 INTERFACE SPLDRS 00085 MODULE PROCEDURE SPLDRS 00086 END INTERFACE 00087 INTERFACE SPLPR 00088 MODULE PROCEDURE SPLPR 00089 END INTERFACE 00090 INTERFACE SPLDV 00091 MODULE PROCEDURE SPLDV 00092 END INTERFACE 00093 INTERFACE SPLD2V 00094 MODULE PROCEDURE SPLD2V 00095 END INTERFACE 00096 INTERFACE SPLPV 00097 MODULE PROCEDURE SPLPV 00098 END INTERFACE 00099 INTERFACE SPLP 00100 MODULE PROCEDURE SPLP 00101 END INTERFACE 00102 INTERFACE SPLTT 00103 MODULE PROCEDURE SPLTT 00104 END INTERFACE 00105 INTERFACE SPLR 00106 MODULE PROCEDURE SPLR_1 00107 MODULE PROCEDURE SPLR_2 00108 END INTERFACE 00109 INTERFACE SPLU 00110 MODULE PROCEDURE SPLU 00111 END INTERFACE 00112 INTERFACE SPLW 00113 MODULE PROCEDURE SPLW 00114 END INTERFACE 00115 INTERFACE SPLVPQ 00116 MODULE PROCEDURE SPLVPQ 00117 END INTERFACE 00118 INTERFACE SPLC 00119 MODULE PROCEDURE SPLC 00120 END INTERFACE 00121 INTERFACE SPLM 00122 MODULE PROCEDURE SPLM 00123 END INTERFACE 00124 INTERFACE SP0NOP 00125 MODULE PROCEDURE SP0NOP 00126 END INTERFACE 00127 INTERFACE SP0CVQ 00128 MODULE PROCEDURE SP0CVQ 00129 END INTERFACE 00130 INTERFACE SPLBSD 00131 MODULE PROCEDURE SPLBSD 00132 END INTERFACE 00133 INTERFACE SPLBSEL 00134 MODULE PROCEDURE SPLBSEL 00135 END INTERFACE 00136 INTERFACE SPLB2C 00137 MODULE PROCEDURE SPLB2C 00138 END INTERFACE 00139 00140 CONTAINS 00141 00142 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00143 00144 !COPIE D'UNE MATRICE DANS UNE AUTRE 00145 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00146 !1. copie de matrices potentiellement de tailles différentes 00147 !SPLR 00148 SUBROUTINE SSCOPY_1(A,B,IA,IB) 00149 ! 00150 IMPLICIT NONE 00151 00152 REAL, DIMENSION(:), INTENT(IN) :: A 00153 INTEGER, INTENT(IN) :: IA 00154 INTEGER, INTENT(IN) :: IB 00155 REAL, DIMENSION(:), INTENT(OUT) :: B 00156 00157 INTEGER :: I, J, K 00158 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00159 00160 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SSCOPY_1',0,ZHOOK_HANDLE) 00161 J = IB 00162 K = IA 00163 DO I=1,SIZE(A) 00164 B(J) = A(K) 00165 J = J + IB 00166 K = K + IA 00167 ENDDO 00168 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SSCOPY_1',1,ZHOOK_HANDLE) 00169 00170 END SUBROUTINE SSCOPY_1 00171 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00172 !1. copie de matrices potentiellement de tailles différentes 00173 SUBROUTINE SSCOPY_2(A,B,IA,IB) 00174 ! 00175 IMPLICIT NONE 00176 00177 REAL, DIMENSION(:,:), INTENT(IN) :: A 00178 INTEGER, INTENT(IN) :: IA 00179 INTEGER, INTENT(IN) :: IB 00180 REAL, DIMENSION(:,:), INTENT(OUT) :: B 00181 00182 REAL,DIMENSION(SIZE(A,1)*SIZE(A,2)) :: A1 00183 REAL,DIMENSION(SIZE(B,1)*SIZE(B,2)) :: B1 00184 00185 INTEGER :: N 00186 INTEGER :: I, J, K 00187 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00188 00189 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SSCOPY_2',0,ZHOOK_HANDLE) 00190 N=SIZE(A,1)*SIZE(A,2) 00191 00192 DO J=1,SIZE(A,2) 00193 A1((J-1)*SIZE(A,1)+1:J*SIZE(A,1))=A(:,J) 00194 ENDDO 00195 00196 J = IB 00197 K = IA 00198 DO I=1,N-1 00199 B1(J) = A1(K) 00200 J = J + IB 00201 K = K + IA 00202 ENDDO 00203 00204 DO J=1,SIZE(B,2) 00205 B(:,J)=B1((J-1)*SIZE(B,1)+1:J*SIZE(B,1)) 00206 ENDDO 00207 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SSCOPY_2',1,ZHOOK_HANDLE) 00208 00209 END SUBROUTINE SSCOPY_2 00210 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00211 00212 !OPERATIONS SUR LES MATRICES 00213 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00214 !2. A(I) B(I*J) R(I*J) R(:,K)=B(:,K)/A(:) 00215 SUBROUTINE MXIDML(A,B,R) 00216 !SPLTT, SPLR 00217 IMPLICIT NONE 00218 00219 REAL, DIMENSION(:), INTENT(IN) :: A 00220 REAL, DIMENSION(:,:), INTENT(IN) :: B 00221 REAL, DIMENSION(:,:), INTENT(OUT):: R 00222 00223 INTEGER :: K 00224 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00225 00226 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MXIDML',0,ZHOOK_HANDLE) 00227 IF (SIZE(B,1).EQ.0 .OR. SIZE(B,2).EQ.0 .AND. LHOOK) & 00228 CALL DR_HOOK('MODE_SPLINES:MXIDML',1,ZHOOK_HANDLE) 00229 IF (SIZE(B,1).EQ.0 .OR. SIZE(B,2).EQ.0) RETURN 00230 00231 DO K=1,SIZE(B,2) 00232 R(:,K)=B(:,K)/A(:) 00233 ENDDO 00234 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MXIDML',1,ZHOOK_HANDLE) 00235 00236 END SUBROUTINE MXIDML 00237 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00238 !3. produit bizarre de matrices 00239 SUBROUTINE MTXAXM(A,B,R) 00240 !SPLR,SPLU,SPLVPQ 00241 IMPLICIT NONE 00242 00243 REAL, DIMENSION(:,:),INTENT(IN) :: A 00244 REAL, DIMENSION(:,:),INTENT(IN) :: B 00245 REAL, DIMENSION(:,:),INTENT(OUT)::R 00246 00247 INTEGER :: I,J,K, L 00248 REAL :: RIJ 00249 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00250 00251 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MTXAXM',0,ZHOOK_HANDLE) 00252 R(:,:)=0. 00253 DO I=1,SIZE(A,2) 00254 DO J=1,SIZE(A,1) 00255 RIJ=0. 00256 DO K=1,SIZE(A,1) 00257 RIJ=RIJ+B(J,K)*A(K,I) 00258 ENDDO 00259 DO L=1,I 00260 R(L,I)=R(L,I)+RIJ*A(J,L) 00261 R(I,L)=R(L,I) 00262 ENDDO 00263 ENDDO 00264 ENDDO 00265 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MTXAXM',1,ZHOOK_HANDLE) 00266 00267 END SUBROUTINE MTXAXM 00268 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00269 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00270 !4. produit simple de matrices (mettre les deux dans la même interface) 00271 SUBROUTINE MXMSPL_1(A,B,R) 00272 !SPLE, SPLC, SPLTT, SPLR, SPLW, SPLVPQ 00273 IMPLICIT NONE 00274 00275 REAL, DIMENSION(:,:),INTENT(IN) :: A 00276 REAL, DIMENSION(:,:),INTENT(IN) :: B 00277 REAL, DIMENSION(:,:),INTENT(OUT)::R 00278 00279 INTEGER :: I,J,K 00280 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00281 00282 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MXMSPL_1',0,ZHOOK_HANDLE) 00283 R(:,:)=0. 00284 00285 DO I=1,SIZE(A,1) 00286 DO J=1,SIZE(B,2) 00287 DO K=1,SIZE(A,2) 00288 R(I,J)=R(I,J)+A(I,K)*B(K,J) 00289 ENDDO 00290 ENDDO 00291 ENDDO 00292 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MXMSPL_1',1,ZHOOK_HANDLE) 00293 00294 END SUBROUTINE MXMSPL_1 00295 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00296 !4. produit simple de matrices 2 00297 SUBROUTINE MXMSPL_2(A,B,R) 00298 ! 00299 IMPLICIT NONE 00300 00301 REAL, DIMENSION(:,:),INTENT(IN) :: A 00302 REAL, DIMENSION(:),INTENT(IN) :: B 00303 REAL, DIMENSION(:),INTENT(OUT)::R 00304 00305 INTEGER :: I,J 00306 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00307 00308 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MXMSPL_2',0,ZHOOK_HANDLE) 00309 R(:)=0. 00310 00311 DO I=1,SIZE(A,1) 00312 DO J=1,SIZE(A,2) 00313 R(I)=R(I)+A(I,J)*B(J) 00314 ENDDO 00315 ENDDO 00316 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MXMSPL_2',1,ZHOOK_HANDLE) 00317 00318 END SUBROUTINE MXMSPL_2 00319 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00320 !4. produit simple de matrices 3 00321 SUBROUTINE MXMSPL_3(A,B,RES) 00322 ! 00323 IMPLICIT NONE 00324 00325 REAL, DIMENSION(:),INTENT(IN) :: A 00326 REAL, DIMENSION(:),INTENT(IN) :: B 00327 REAL,INTENT(OUT)::RES 00328 00329 INTEGER :: I 00330 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00331 00332 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MXMSPL_3',0,ZHOOK_HANDLE) 00333 RES=0. 00334 00335 DO I=1,SIZE(A) 00336 RES=RES+A(I)*B(I) 00337 ENDDO 00338 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:MXMSPL_3',1,ZHOOK_HANDLE) 00339 00340 END SUBROUTINE MXMSPL_3 00341 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00342 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00343 !5. inversion de la matrice A. 00344 SUBROUTINE SMXINV(A,IREP) 00345 !SPLTT 00346 IMPLICIT NONE 00347 00348 REAL, DIMENSION(:,:), INTENT(INOUT) :: A 00349 INTEGER, INTENT(OUT) :: IREP 00350 00351 REAL, DIMENSION(SIZE(A,1)*SIZE(A,2)) :: A2 00352 INTEGER, DIMENSION(SIZE(A,1)) :: INDEX 00353 REAL, DIMENSION(SIZE(A,1)) :: RI 00354 INTEGER :: N, NP, NP1, I, J, JJ, K, L, KK, KJ, JL 00355 INTEGER :: IJ0, JI0, IJ, JI 00356 REAL :: PIVOT, ELM 00357 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00358 00359 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SMXINV',0,ZHOOK_HANDLE) 00360 IREP=0 00361 N=SIZE(A,1) 00362 NP1=N+1 00363 00364 INDEX(:)=1 00365 00366 DO J=1,N 00367 !les colonnes sont à la suite 00368 A2((J-1)*N+1:J*N)=A(:,J) 00369 ENDDO 00370 00371 DO I=1,N 00372 !find pivot 00373 PIVOT=0. 00374 JJ=1 00375 !dans la diagonale, on retient le plus grand élément 00376 DO J=1,N 00377 IF(INDEX(J).NE.0) THEN 00378 ELM=ABS(A2(JJ)) 00379 IF (ELM.GT.PIVOT) THEN 00380 PIVOT=ELM 00381 K=J !numéro de la ligne / colonne 00382 !numéro de l'élément dans A2 00383 KK=JJ 00384 ENDIF 00385 ENDIF 00386 JJ=JJ+NP1 00387 ENDDO 00388 INDEX(K)=0 00389 IF (PIVOT.EQ.0.) THEN 00390 IREP=-1 00391 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SMXINV',1,ZHOOK_HANDLE) 00392 RETURN 00393 ENDIF 00394 PIVOT=-A2(KK) 00395 !elimination 00396 KJ=K 00397 NP=N 00398 DO J=1,N 00399 !KJ=K, KJ=K+N, KJ=K+2N ... : avant l'élément pivot, 00400 !on se déplace sur les colonnes 00401 !J=K=3: KJ=K+2N+1, KJ=K+2N+2....!une fois qu'on a passé 00402 !le pivot, on se déplace sur la ligne, à droite de l'élément pivot 00403 !si on est sur le premier élément de la colonne pivot 00404 IF (J==K) THEN 00405 !KJ réfère à l'élément diagonal du pivot 00406 A2(KJ)=1./PIVOT 00407 RI(J)=0. 00408 NP=1 00409 ELSE 00410 !si on est sur un autre élément de du pivot 00411 ELM=-A2(KJ) 00412 RI(J)=ELM/PIVOT 00413 IF (ELM.NE.0.) THEN 00414 JL=J 00415 !JL prend pour valeur les élements de la colonne 00416 !du premier en haut au diagonal -1 00417 DO L=1,J 00418 !à chaque fois A2 de cet élément est incrémenté 00419 A2(JL)=A2(JL)+ELM*RI(L) 00420 JL=JL+N 00421 ENDDO 00422 ENDIF 00423 !élement véritablement en cours 00424 A2(KJ)=RI(J) 00425 ENDIF 00426 !KJ=K+JN si on n'a pas encore passé l'élément de la colonne pivot 00427 !KJ=K+J0N+J si on l'a passé 00428 KJ=KJ+NP 00429 ENDDO 00430 ENDDO 00431 00432 !change the sign and provisional fill-up 00433 DO J=1,N 00434 A(:,J)= A2((J-1)*N+1:J*N) 00435 ENDDO 00436 00437 DO I=1,N 00438 DO J=1,I 00439 A(I,J)=-A(I,J) 00440 A(J,I)=A(I,J) 00441 ENDDO 00442 ENDDO 00443 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SMXINV',1,ZHOOK_HANDLE) 00444 00445 END SUBROUTINE SMXINV 00446 00447 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00448 ! CALCULS PLUS SPECIFIQUES 00449 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00450 !6. IE(ND,M): calcul de la matrice IE en fonction de NDEG. 00451 SUBROUTINE SPLIE(NDEG,IE) 00452 !SPLE, SPLR 00453 IMPLICIT NONE 00454 00455 INTEGER, INTENT(IN) :: NDEG 00456 INTEGER, DIMENSION(:,:), INTENT(OUT) :: IE 00457 00458 INTEGER, DIMENSION(SIZE(IE,1)) :: NV 00459 INTEGER :: I, J, K, T, N, N0, NC 00460 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00461 00462 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLIE',0,ZHOOK_HANDLE) 00463 NV(:)=1 00464 IE(:,1)=0 00465 N=1 00466 N0=1 00467 00468 !boucle d'itération 00469 DO T=1,NDEG 00470 !boucle sur les indices 00471 DO I=1,SIZE(IE,1) 00472 NC=0 00473 !NV(I) varie d'itération T en itération 00474 DO K=NV(I),N0 00475 N=N+1 00476 NC=NC+1 00477 !boucle sur les indices 00478 DO J=1,SIZE(IE,1) 00479 IE(J,N)=IE(J,K) 00480 IF (J.EQ.I) IE(J,N)=IE(J,N)+1 00481 ENDDO 00482 ENDDO 00483 NV(I)=N-NC+1 00484 ENDDO 00485 N0=N 00486 ENDDO 00487 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLIE',1,ZHOOK_HANDLE) 00488 00489 END SUBROUTINE SPLIE 00490 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00491 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00492 !7. calcul de RK en fonction des autres RK(N1,N2) 00493 SUBROUTINE SPLK_1(NORD,X1,X2,G,RK) 00494 !SPLE,SPLU, 00495 IMPLICIT NONE 00496 00497 INTEGER, INTENT(IN) :: NORD 00498 REAL, DIMENSION(:), INTENT(IN) :: X1 00499 REAL, DIMENSION(:,:), INTENT(IN) :: X2 00500 REAL, DIMENSION(:,:), INTENT(IN) :: G 00501 REAL, DIMENSION(:), INTENT(OUT):: RK 00502 00503 INTEGER :: I, ID, JD, EXP, ISIGNE 00504 INTEGER :: ND, N, DI, DJ 00505 REAL :: D2 00506 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00507 00508 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLK_1',0,ZHOOK_HANDLE) 00509 ND=SIZE(X1) 00510 N=SIZE(X2,2) 00511 00512 !Exposant du noyau 00513 EXP=(2*NORD-ND)/2. 00514 00515 IF (MOD(ND,2).EQ.0) THEN 00516 !Calcul de K dans le cas d'un espace de dimension paire 00517 ISIGNE=(-1)**(1+NORD+ND/2) 00518 ELSE 00519 !Calcul de K dans le cas d'un espace de dimension impair 00520 ISIGNE=(-1)**(NORD+ND/2) 00521 ENDIF 00522 00523 RK(:)=0. 00524 00525 DO I=1,N 00526 D2=0. 00527 DO ID=1,ND 00528 DO JD=1,ND 00529 DI=X1(ID)-X2(ID,I) 00530 DJ=X1(JD)-X2(JD,I) 00531 D2=D2+G(ID,JD)*DI*DJ 00532 ENDDO 00533 ENDDO 00534 IF (D2.NE.0.) THEN 00535 RK(I) = ISIGNE*D2**EXP 00536 IF (MOD(ND,2).EQ.0) RK(I)=RK(I)*0.5*LOG(D2) 00537 ENDIF 00538 ENDDO 00539 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLK_1',1,ZHOOK_HANDLE) 00540 00541 END SUBROUTINE SPLK_1 00542 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00543 !7. calcul de RK en fonction des autres RK(N1,N2) 00544 SUBROUTINE SPLK_2(NORD,X1,X2,G,RK) 00545 00546 IMPLICIT NONE 00547 00548 INTEGER, INTENT(IN) :: NORD 00549 REAL, DIMENSION(:,:), INTENT(IN) :: X1 00550 REAL, DIMENSION(:,:), INTENT(IN) :: X2 00551 REAL, DIMENSION(:,:), INTENT(IN) :: G 00552 REAL, DIMENSION(:,:), INTENT(OUT):: RK 00553 00554 INTEGER :: I1, I2, ID, JD, ISIGNE 00555 INTEGER :: ND, N1, N2 00556 REAL :: EXP, D2, DI, DJ 00557 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00558 00559 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLK_2',0,ZHOOK_HANDLE) 00560 ND=SIZE(X1,1) 00561 N1=SIZE(X1,2) 00562 N2=SIZE(X2,2) 00563 00564 !Exposant du noyau 00565 EXP=(2*NORD-ND)/2. 00566 00567 IF (MOD(ND,2).EQ.0) THEN 00568 !Calcul de K dans le cas d'un espace de dimension paire 00569 ISIGNE=(-1)**(1+NORD+ND/2) 00570 ELSE 00571 !Calcul de K dans le cas d'un espace de dimension impair 00572 ISIGNE=(-1)**(NORD+ND/2) 00573 ENDIF 00574 00575 RK(:,:)=0. 00576 00577 DO I1=1,N1 00578 DO I2=1,N2 00579 D2=0. 00580 DO ID=1,ND 00581 DO JD=1,ND 00582 DI=X1(ID,I1)-X2(ID,I2) 00583 DJ=X1(JD,I1)-X2(JD,I2) 00584 D2=D2+G(ID,JD)*DI*DJ 00585 ENDDO 00586 ENDDO 00587 IF (D2.NE.0.) THEN 00588 RK(I1,I2) = ISIGNE*D2**EXP 00589 IF (MOD(ND,2).EQ.0) RK(I1,I2)=RK(I1,I2)*0.5*LOG(D2) 00590 ENDIF 00591 ENDDO 00592 ENDDO 00593 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLK_2',1,ZHOOK_HANDLE) 00594 00595 END SUBROUTINE SPLK_2 00596 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00597 !8. on cherche l'indice I de la 3ème dimension de AI tel que XE(ID) soit compris 00598 !entre AI(ID,1,I) et AI(ID,2,I+1). Il est compris entre 1 et NSD(ID), 00599 !par défaut il est égal à NDSD(ID) 00600 SUBROUTINE SPLBFIN(ID,NSD,XE,AI,IDX,IDX1) 00601 !SPLB2E 00602 IMPLICIT NONE 00603 00604 INTEGER, INTENT(IN) :: ID 00605 INTEGER, DIMENSION(:), INTENT(IN) :: NSD 00606 REAL, DIMENSION(:), INTENT(IN) :: XE 00607 REAL, DIMENSION(:,:,:), INTENT(IN) :: AI 00608 INTEGER, INTENT(OUT) :: IDX 00609 INTEGER, INTENT(OUT) :: IDX1 00610 00611 !INTEGER, PARAMETER: NSDMAX=20 00612 INTEGER :: I, I1 00613 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00614 00615 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBFIN',0,ZHOOK_HANDLE) 00616 IDX1=0 00617 DO I=1,NSD(ID)-1 00618 IF (XE(ID).LE.AI(ID,2,I)) THEN 00619 IDX=I 00620 I1=I+1 00621 IF (XE(ID).GT.AI(ID,1,I1)) IDX1=I1 00622 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBFIN',1,ZHOOK_HANDLE) 00623 RETURN 00624 ENDIF 00625 ENDDO 00626 IDX=NSD(ID) 00627 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBFIN',1,ZHOOK_HANDLE) 00628 00629 END SUBROUTINE SPLBFIN 00630 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00631 !9. X(ND,N) IE(ND,M) : T(N,M) T=produit(X**IE) quand IE est non nul 00632 SUBROUTINE SPLT_1(X,IE,T) 00633 !SPLE, SPLTT 00634 IMPLICIT NONE 00635 00636 REAL, DIMENSION(:), INTENT(IN) :: X 00637 INTEGER, DIMENSION(:,:), INTENT(IN) :: IE 00638 REAL, DIMENSION(:), INTENT(OUT) :: T 00639 ! 00640 INTEGER :: J,K 00641 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00642 00643 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLT_1',0,ZHOOK_HANDLE) 00644 T(:)=1. 00645 00646 DO J=1,SIZE(IE,2) 00647 DO K=1,SIZE(X,1) 00648 IF (IE(K,J).NE.0) T(J)=T(J)*X(K)**IE(K,J) 00649 ENDDO 00650 ENDDO 00651 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLT_1',1,ZHOOK_HANDLE) 00652 00653 END SUBROUTINE SPLT_1 00654 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00655 !9. X(ND,N) IE(ND,M) : T(N,M) T=produit(X**IE) quand IE est non nul 00656 SUBROUTINE SPLT_2(X,IE,T) 00657 ! 00658 IMPLICIT NONE 00659 00660 REAL, DIMENSION(:,:), INTENT(IN) :: X 00661 INTEGER, DIMENSION(:,:), INTENT(IN) :: IE 00662 REAL, DIMENSION(:,:), INTENT(OUT) :: T 00663 ! 00664 INTEGER :: I,J,K 00665 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00666 00667 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLT_2',0,ZHOOK_HANDLE) 00668 T(:,:)=1. 00669 00670 DO J=1,SIZE(IE,2) 00671 DO I=1,SIZE(X,2) 00672 DO K=1,SIZE(X,1) 00673 IF (IE(K,J).NE.0) T(I,J)=T(I,J)*X(K,I)**IE(K,J) 00674 ENDDO 00675 ENDDO 00676 ENDDO 00677 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLT_2',1,ZHOOK_HANDLE) 00678 00679 END SUBROUTINE SPLT_2 00680 00681 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00682 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00683 ! PARTIE "E" 00684 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00685 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00686 !10. 00687 SUBROUTINE SPLE(NORD,X,G,C,XE,ZE) 00688 !SPLB2E1 00689 IMPLICIT NONE 00690 00691 INTEGER, INTENT(IN) :: NORD 00692 REAL, DIMENSION(:,:),INTENT(IN) :: X !ND,N 00693 REAL, DIMENSION(:,:), INTENT(IN) :: G !ND,ND 00694 REAL, DIMENSION(:),INTENT(IN) :: C ! N+M 00695 REAL, DIMENSION(:),INTENT(IN) :: XE !ND 00696 REAL, INTENT(OUT) :: ZE 00697 00698 INTEGER, DIMENSION(SIZE(X,1),SIZE(C)-SIZE(X,2)) :: IW !ND,M 00699 REAL, DIMENSION(SIZE(X,2)) :: WE1 00700 REAL, DIMENSION(SIZE(IW,2)) :: WE2 00701 REAL :: ZOUT1,ZOUT2 00702 INTEGER :: M, N, NE, M1, M2, M3 00703 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00704 00705 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLE',0,ZHOOK_HANDLE) 00706 N=SIZE(X,2) 00707 M=SIZE(C)-N 00708 00709 ! Calcul de la matrice Kxe 00710 CALL SPLK(NORD,XE,X,G,WE1) 00711 ! Calcul de Kxe.c 00712 CALL MXMSPL(WE1,C(1:N),ZOUT1) 00713 00714 IF (M.NE.0) THEN 00715 ! Generation des exposants des monomes 00716 CALL SPLIE(NORD-1,IW) !IW(ND,M) 00717 ! Calcul de la matrice Txe 00718 CALL SPLT(XE,IW,WE2) 00719 ! Calcul Txe.d 00720 CALL MXMSPL(WE2,C(N+1:N+M),ZOUT2) 00721 ! Calcul de Txe.d + Kxe.c 00722 ZE = ZOUT1+ZOUT2 00723 ENDIF 00724 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLE',1,ZHOOK_HANDLE) 00725 00726 END SUBROUTINE SPLE 00727 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00728 !11. 00729 SUBROUTINE SPLB2E1(NORD,M,G,C,XE,IDX,IDX1,IDY,IDY1,ZE) 00730 !SPLB2E 00731 00732 IMPLICIT NONE 00733 00734 INTEGER, INTENT(IN) :: NORD 00735 INTEGER, INTENT(IN) :: M 00736 REAL, DIMENSION(:,:),INTENT(IN) :: G 00737 REAL, DIMENSION(:,:,:), INTENT(IN) :: C 00738 REAL, DIMENSION(:), INTENT(IN) :: XE 00739 INTEGER, INTENT(IN) :: IDX,IDX1,IDY,IDY1 00740 REAL, INTENT(OUT) :: ZE 00741 00742 INTEGER :: NN,NM 00743 REAL :: RX,RY,ALPHAX,ALPHAY,SUM_AX,SUM_AY,SUM_A 00744 REAL :: ZXY, ZX1Y, ZXY1, ZX1Y1 00745 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00746 00747 00748 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLB2E1',0,ZHOOK_HANDLE) 00749 00750 RX=0 00751 RY=0 00752 IF (IDX1.NE.0) RX=(XE(1)-AI(1,1,IDX1))/(AI(1,2,IDX)-AI(1,1,IDX1)) 00753 IF (IDY1.NE.0) RY=(XE(2)-AI(2,1,IDY1))/(AI(2,2,IDY)-AI(2,1,IDY1)) 00754 ALPHAX=(RX-1.)*(RX-1.)*(2*RX+1.) 00755 ALPHAY=(RY-1.)*(RY-1.)*(2*RY+1.) 00756 SUM_AX=0. 00757 SUM_AY=0. 00758 SUM_A=0. 00759 ZXY=0. 00760 ZX1Y=0. 00761 ZXY1=0. 00762 ZX1Y1=0. 00763 00764 NM=0 00765 00766 IF (NS(IDX,IDY).GE.M) THEN 00767 NM=1 00768 SUM_AX=SUM_AX+ALPHAX 00769 SUM_AY=SUM_AY+ALPHAY 00770 SUM_A=SUM_A+ALPHAX*ALPHAY 00771 NN=NS(IDX,IDY) 00772 !N=NN 00773 CALL SPLE(NORD,XS(:,1:NN,IDX,IDY),G,C(1:NN+M,IDX,IDY),XE,ZXY) 00774 ENDIF 00775 00776 IF (IDX1.NE.0) THEN 00777 IF (NS(IDX1,IDY).GE.M) THEN 00778 NM=1 00779 SUM_AX=SUM_AX+1.-ALPHAX 00780 SUM_A=SUM_A+(1-ALPHAX)*ALPHAY 00781 NN=NS(IDX1,IDY) 00782 CALL SPLE(NORD,XS(:,1:NN,IDX1,IDY),G,C(1:NN+M,IDX1,IDY),XE,ZX1Y) 00783 ENDIF 00784 ENDIF 00785 00786 IF (IDY1.NE.0) THEN 00787 IF (NS(IDX,IDY1).GE.M) THEN 00788 NM=1 00789 SUM_AY=SUM_AY+1.-ALPHAY 00790 SUM_A=SUM_A+ALPHAX*(1-ALPHAY) 00791 NN=NS(IDX,IDY1) 00792 CALL SPLE(NORD,XS(:,1:NN,IDX,IDY1),G,C(1:NN+M,IDX,IDY1),XE,ZXY1) 00793 ENDIF 00794 ENDIF 00795 00796 IF (IDX1.NE.0 .AND. IDY1.NE.0) THEN 00797 IF (NS(IDX1,IDY1).GE.M) THEN 00798 NM=1 00799 SUM_A=SUM_A+(1-ALPHAX)*(1-ALPHAY) 00800 NN=NS(IDX1,IDY1) 00801 CALL SPLE(NORD,XS(:,1:NN,IDX1,IDY1),G,C(1:NN+M,IDX1,IDY1),XE,ZX1Y1) 00802 ENDIF 00803 ENDIF 00804 00805 IF (NM==1) THEN 00806 IF (IDX1.NE.0 .AND. IDY1.EQ.0) THEN 00807 ZE=(ALPHAX*ZXY+(1-ALPHAX)*ZX1Y)/SUM_AX 00808 ELSEIF (IDX1.EQ.0 .AND. IDY1.NE.0) THEN 00809 ZE=(ALPHAY*ZXY+(1-ALPHAY)*ZXY1)/SUM_AY 00810 ELSEIF (IDX1.NE.0 .AND. IDY1.NE.0) THEN 00811 ZE=(ALPHAX*ALPHAY *ZXY + (1-ALPHAX)* ALPHAY*ZX1Y + & 00812 ALPHAX*(1-ALPHAY)*ZXY1 + (1-ALPHAX)*(1-ALPHAY)*ZX1Y1) / SUM_A 00813 ELSE 00814 ZE=ZXY 00815 ENDIF 00816 ENDIF 00817 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLB2E1',1,ZHOOK_HANDLE) 00818 00819 END SUBROUTINE SPLB2E1 00820 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00821 ! APPELEE PAR INTERPOL_SPLINES.F90 00822 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00823 !12. 00824 SUBROUTINE SPLB2E(NORD,M,NSDI,NSDJ,G,C,XE,ZE) 00825 00826 IMPLICIT NONE 00827 00828 INTEGER, INTENT(IN) :: NORD 00829 INTEGER, INTENT(IN) :: M 00830 INTEGER, INTENT(IN) :: NSDI,NSDJ 00831 REAL, DIMENSION(:,:), INTENT(IN) :: G 00832 REAL,DIMENSION(:,:,:), INTENT(IN) :: C 00833 REAL, DIMENSION(:,:), INTENT(IN) :: XE 00834 REAL, DIMENSION(:), INTENT(OUT) :: ZE 00835 00836 INTEGER, DIMENSION(2) :: NSD 00837 INTEGER :: J, ID, NE 00838 INTEGER :: IDX, IDX1, IDY, IDY1 00839 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00840 00841 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLB2E',0,ZHOOK_HANDLE) 00842 NE=SIZE(XE,2) 00843 00844 NSD(1)=NSDI 00845 NSD(2)=NSDJ 00846 00847 DO J=1,NE 00848 ID=1 00849 CALL SPLBFIN(ID,NSD,XE(:,J),AI,IDX,IDX1) 00850 ID=2 00851 CALL SPLBFIN(ID,NSD,XE(:,J),AI,IDY,IDY1) 00852 CALL SPLB2E1(NORD,M,G,C,XE(:,J),IDX,IDX1,IDY,IDY1,ZE(J)) 00853 ENDDO 00854 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLB2E',1,ZHOOK_HANDLE) 00855 00856 END SUBROUTINE SPLB2E 00857 00858 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00859 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00860 ! PARTIE "C" 00861 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00862 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00863 00864 ! Calculs EISRS1 00865 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00866 !13. 00867 SUBROUTINE TRED2(A,D,E,Z) 00868 !EISRS1 00869 IMPLICIT NONE 00870 00871 REAL, DIMENSION(:,:), INTENT(IN) :: A 00872 REAL, DIMENSION(:), INTENT(OUT) :: D 00873 REAL, DIMENSION(:), INTENT(OUT) :: E 00874 REAL, DIMENSION(:,:), INTENT(OUT):: Z 00875 00876 INTEGER :: N, I, II, J, K, L, JP1 00877 REAL :: F, G, H, SCALE, HH 00878 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00879 00880 00881 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:TRED2',0,ZHOOK_HANDLE) 00882 N=SIZE(A,1) 00883 00884 !Z=A sur la partie gauceh jusqu'à la diagonale 00885 DO I=1,N 00886 DO J=1,I 00887 Z(I,J)=A(I,J) 00888 ENDDO 00889 ENDDO 00890 00891 IF (N.NE.1) THEN 00892 DO II=2,N 00893 !I varie de N à 2 00894 I=N+2-II 00895 !L varie de N-1 à 1 00896 L=I-1 00897 H=0.0 00898 SCALE=0.0 00899 IF (L.GE.2) THEN 00900 !scale = somme des éléments de la ligne à l'exclusion du diagonal 00901 DO K=1,L 00902 SCALE=SCALE+ABS(Z(I,K)) 00903 ENDDO 00904 IF (SCALE.NE.0.0) THEN 00905 !normalisation de Z par la valeur de la ligne 00906 DO K=1,L 00907 Z(I,K)=Z(I,K)/SCALE 00908 !H = somme des carrés de Z / (somme des ABS(Z)) au carré 00909 H=H+Z(I,K)**2 00910 ENDDO 00911 F=Z(I,L) 00912 !SQRT(H) signé par lélément diagonal-1 de la ligne 00913 G=-SIGN(SQRT(H),F) 00914 !moyenne des éléments au carrés de la ligne, signée 00915 E(I)=SCALE*G 00916 H=H-F*G 00917 Z(I,L)=F-G 00918 F=0.0 00919 DO J=1,L 00920 Z(J,I)=Z(I,J)/(SCALE*H) 00921 G=0.0 00922 DO K=1,J 00923 G=G+Z(J,K)*Z(I,K) 00924 ENDDO 00925 JP1=J+1 00926 IF (L.GE.JP1) THEN 00927 DO K=JP1,L 00928 G=G+Z(K,J)*Z(I,K) 00929 ENDDO 00930 ENDIF 00931 E(J)=G/H 00932 F=F+E(J)*Z(I,J) 00933 ENDDO 00934 HH=F/(H+H) 00935 DO J=1,L 00936 F=Z(I,J) 00937 G=E(J)-HH*F 00938 E(J)=G 00939 DO K=1,J 00940 Z(J,K)=Z(J,K)-F*E(K)-G*Z(I,K) 00941 ENDDO 00942 ENDDO 00943 DO K=1,L 00944 Z(I,K)=SCALE*Z(I,K) 00945 ENDDO 00946 ELSE 00947 E(I)=Z(I,L) 00948 ENDIF 00949 ELSE 00950 E(I)=Z(I,L) 00951 ENDIF 00952 D(I)=H 00953 ENDDO 00954 ENDIF 00955 00956 D(1)=0.0 00957 E(1)=0.0 00958 DO I=1,N 00959 L=I-1 00960 IF (D(I).NE.0.0) THEN 00961 DO J=1,L 00962 G=0.0 00963 DO K=1,L 00964 G=G+Z(I,K)*Z(K,J) 00965 ENDDO 00966 DO K=1,L 00967 Z(K,J)=Z(K,J)-G*Z(K,I) 00968 ENDDO 00969 ENDDO 00970 ENDIF 00971 D(I)=Z(I,I) 00972 Z(I,I)=1.0 00973 IF (L.GE.1) THEN 00974 DO J=1,L 00975 Z(I,J)=0.0 00976 Z(J,I)=0.0 00977 ENDDO 00978 ENDIF 00979 ENDDO 00980 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:TRED2',1,ZHOOK_HANDLE) 00981 00982 END SUBROUTINE TRED2 00983 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 00984 !14. 00985 SUBROUTINE TQL2_2(D,E,Z,IERR) 00986 !EISRS1 00987 IMPLICIT NONE 00988 00989 REAL, DIMENSION(:), INTENT(INOUT) :: D 00990 REAL, DIMENSION(:), INTENT(INOUT) :: E 00991 REAL, DIMENSION(:,:), INTENT(OUT) :: Z 00992 INTEGER, INTENT(OUT) :: IERR 00993 00994 INTEGER :: I, II, J, K, L, M, N, MML 00995 REAL :: B, C, F, G, H, P, R, S 00996 REAL :: MACHEP 00997 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00998 00999 01000 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:TQL2_2',0,ZHOOK_HANDLE) 01001 MACHEP=2.**(-52) 01002 IERR=0 01003 N=SIZE(D) 01004 01005 IF (N.NE.1) THEN 01006 DO I=2,N 01007 E(I-1)=E(I) 01008 ENDDO 01009 F=0.0 01010 B=0.0 01011 E(N)=0.0 01012 DO L=1,N 01013 J=0 01014 H=MACHEP*(ABS(D(L))+ABS(E(L))) 01015 IF (B.LT.H) B=H 01016 DO M=L,N 01017 IF (ABS(E(M)).LE.B) EXIT 01018 ENDDO 01019 IF (M.NE.L) THEN 01020 DO WHILE (ABS(E(L)).GT.B) 01021 IF (J.EQ.30) THEN 01022 IERR=L 01023 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:TQL2_2',1,ZHOOK_HANDLE) 01024 RETURN 01025 ENDIF 01026 J=J+1 01027 P=(D(L+1)-D(L)) / (2.0 * E(L)) 01028 R=SQRT(P*P+1.0) 01029 H=D(L)-E(L) / (P+SIGN(R,P)) 01030 DO I=L,N 01031 D(I)=D(I)-H 01032 ENDDO 01033 F=F+H 01034 P=D(M) 01035 C=1.0 01036 S=0.0 01037 MML=M-L 01038 DO II=1,MML 01039 I=M-II 01040 G=C*E(I) 01041 H=C*P 01042 IF (ABS(P).GE.ABS(E(I))) THEN 01043 C=E(I)/P 01044 R=SQRT(C*C+1.0) 01045 E(I+1)=S*P*R 01046 S=C/R 01047 C=1.0/R 01048 ELSE 01049 C=P/E(I) 01050 R=SQRT(C*C+1.0) 01051 E(I+1)=S*E(I)*R 01052 S=1.0/R 01053 C=C*S 01054 ENDIF 01055 P=C*D(I)-S*G 01056 D(I+1)=H+S*(C*G+S*D(I)) 01057 DO K=1,N 01058 H=Z(K,I+1) 01059 Z(K,I+1)=S*Z(K,I)+C*H 01060 Z(K,I)=C*Z(K,I)-S*H 01061 ENDDO 01062 ENDDO 01063 E(L)=S*P 01064 D(L)=C*P 01065 ENDDO 01066 ENDIF 01067 D(L)=D(L)+F 01068 ENDDO 01069 DO II=2,N 01070 I=II-1 01071 K=I 01072 P=D(I) 01073 DO J=II,N 01074 IF (D(J).LT.P) THEN 01075 K=J 01076 P=D(J) 01077 ENDIF 01078 ENDDO 01079 IF (K.NE.I) THEN 01080 D(K)=D(I) 01081 D(I)=P 01082 DO J=1,N 01083 P=Z(J,I) 01084 Z(J,I)=Z(J,K) 01085 Z(J,K)=P 01086 ENDDO 01087 ENDIF 01088 ENDDO 01089 ENDIF 01090 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:TQL2_2',1,ZHOOK_HANDLE) 01091 01092 END SUBROUTINE TQL2_2 01093 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01094 !15. SPLR,SPLU 01095 SUBROUTINE EISRS1(AR,WR,ZR,WORK,IERR) 01096 ! 01097 ! all eigenvalues and corresponding eigenvectors of a real 01098 ! symmetric matrix 01099 ! 01100 IMPLICIT NONE 01101 !toutes les dimensions sont N: ok 01102 REAL, DIMENSION(:,:), INTENT(IN) :: AR 01103 REAL, DIMENSION(:), INTENT(OUT) :: WR 01104 REAL, DIMENSION(:,:), INTENT(OUT) :: ZR 01105 REAL, DIMENSION(:), INTENT(OUT) :: WORK 01106 INTEGER, INTENT(OUT) :: IERR 01107 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01108 01109 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:EISRS1',0,ZHOOK_HANDLE) 01110 CALL TRED2(AR,WR,WORK,ZR) 01111 CALL TQL2_2(WR,WORK,ZR,IERR) 01112 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:EISRS1',1,ZHOOK_HANDLE) 01113 01114 END SUBROUTINE EISRS1 01115 01116 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01117 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01118 ! GROUPES DES PROCEDURES V 01119 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01120 !16. 01121 SUBROUTINE SPLV(B,W,N,P,RES) 01122 !SP0NOP 01123 IMPLICIT NONE 01124 01125 REAL,DIMENSION(:),INTENT(IN) :: B 01126 REAL,DIMENSION(:), INTENT(IN) :: W 01127 INTEGER, INTENT(IN) :: N 01128 REAL, INTENT(IN) :: P 01129 REAL, INTENT(OUT) :: RES 01130 ! 01131 INTEGER :: A, D, I 01132 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01133 01134 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLV',0,ZHOOK_HANDLE) 01135 A=0. 01136 RES=0. 01137 DO I=1,SIZE(B) 01138 D=B(I)+N*P 01139 A=A+1./D 01140 RES=RES+W(I)**2 / D**2 01141 ENDDO 01142 RES=N*RES/A**2 01143 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLV',1,ZHOOK_HANDLE) 01144 01145 END SUBROUTINE SPLV 01146 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01147 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01148 ! GROUPE DES PROCEDURES P 01149 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01150 !PROCEDURES INTERMEDIAIRES POUR SPLPS2V 01151 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01152 !17. S2VI=somme(W(i)**2)/size(W) 01153 SUBROUTINE SPLS2VI(W,RES) 01154 ! 01155 IMPLICIT NONE 01156 01157 REAL, DIMENSION(:), INTENT(IN) :: W 01158 REAL, INTENT(OUT) :: RES 01159 01160 INTEGER :: I 01161 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01162 01163 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLS2VI',0,ZHOOK_HANDLE) 01164 RES = 0. 01165 DO I=1,SIZE(W) 01166 RES = RES + W(I)**2 01167 ENDDO 01168 RES=RES/SIZE(W) 01169 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLS2VI',1,ZHOOK_HANDLE) 01170 01171 END SUBROUTINE SPLS2VI 01172 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01173 !18. calcul de BN en fonction de W, B, N et P 01174 SUBROUTINE SPLBVM(B,W,N,P,RES) 01175 !SPLS2V, SP0VPQ 01176 IMPLICIT NONE 01177 01178 REAL, DIMENSION(:), INTENT(IN) :: B, W 01179 INTEGER, INTENT(IN) :: N 01180 REAL, INTENT(IN) :: P 01181 REAL, INTENT(OUT) :: RES 01182 ! 01183 REAL :: A 01184 INTEGER :: I 01185 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01186 01187 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBVM',0,ZHOOK_HANDLE) 01188 RES=0. 01189 A=0. 01190 DO I=1,SIZE(B) 01191 RES=RES+(W(I)**2)/(B(I)+N*P)**2 01192 A=A+1./(B(I)+N*P) 01193 ENDDO 01194 RES=RES/A 01195 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBVM',1,ZHOOK_HANDLE) 01196 01197 END SUBROUTINE SPLBVM 01198 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01199 !19. 01200 SUBROUTINE SPLS2V(BI,WI,N,P,RES) 01201 !SPLPS2V, SPLP 01202 IMPLICIT NONE 01203 01204 REAL, DIMENSION(:), INTENT(IN) :: BI !dimension N-M 01205 REAL, DIMENSION(:), INTENT(IN) :: WI !dimension N-M 01206 INTEGER, INTENT(IN) :: N 01207 REAL, INTENT(IN) :: P 01208 REAL, INTENT(OUT) :: RES 01209 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01210 01211 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLS2V',0,ZHOOK_HANDLE) 01212 CALL SPLBVM(BI,WI,N,P,RES) 01213 RES=N*P*RES 01214 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLS2V',1,ZHOOK_HANDLE) 01215 01216 END SUBROUTINE SPLS2V 01217 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01218 !20. 01219 SUBROUTINE SPLDS2V(B,W,N,P,RES) 01220 !SPLPS2V 01221 IMPLICIT NONE 01222 01223 REAL, DIMENSION(:), INTENT(IN) :: B 01224 REAL, DIMENSION(:), INTENT(IN) :: W 01225 INTEGER, INTENT(IN) :: N 01226 REAL, INTENT(IN) :: P 01227 REAL, INTENT(OUT) :: RES 01228 01229 REAL :: DA, DB, A1, A2, A3, B2, C1, C2 01230 INTEGER :: J, I 01231 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01232 01233 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLDS2V',0,ZHOOK_HANDLE) 01234 DA=0. 01235 DB=0. 01236 DO I=1,SIZE(B) 01237 A1=1./(B(I)+N*P) 01238 A2=A1**2 01239 A3=A1**3 01240 DA=DA+A1 01241 B2=W(J)**2 01242 DO J=1,SIZE(B) 01243 C1=1./(B(J)+N*P) 01244 C2=C1**2 01245 DB=DB+B2*(A2*C1+N*P*(-2*A3*C1+A2*C2)) 01246 ENDDO 01247 ENDDO 01248 RES=(N/(DA**2))*DB 01249 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLDS2V',1,ZHOOK_HANDLE) 01250 01251 END SUBROUTINE SPLDS2V 01252 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01253 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01254 !21. SPLPS2V procedure 01255 SUBROUTINE SPLPS2V(B,W,N,P0,S2,P,IREP) 01256 !SPLP 01257 IMPLICIT NONE 01258 01259 REAL, DIMENSION(:), INTENT(IN) :: B 01260 REAL, DIMENSION(:), INTENT(IN) :: W 01261 INTEGER, INTENT(IN) :: N 01262 REAL, INTENT(INOUT) :: P0 01263 REAL, INTENT(IN) :: S2 01264 REAL, INTENT(OUT) :: P 01265 INTEGER, INTENT(OUT) :: IREP 01266 01267 INTEGER, PARAMETER :: EPS=1.E-2 01268 REAL :: P1, RINF, RS, DRS 01269 INTEGER :: NITER 01270 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01271 01272 ! Calcul de la valeur du s2 de Wahba pour p=infini 01273 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPS2V',0,ZHOOK_HANDLE) 01274 CALL SPLS2VI(W,RINF) 01275 IF (RINF.LE.S2) THEN 01276 IREP=-6 01277 P0=0. 01278 WRITE(*,FMT='(A53,G15.5)') & 01279 "SPLPS2V : S2 > VALEUR DE LA FONCTION POUR P INFINI = ",RINF 01280 P=(B(SIZE(B))**2)/B(1) 01281 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPS2V',1,ZHOOK_HANDLE) 01282 RETURN 01283 ENDIF 01284 01285 P1=0. 01286 ! Iterations de Newton 01287 DO NITER=1,NMAX 01288 P0=P1 01289 CALL SPLS2V(B,W,N,P0,RS) 01290 CALL SPLDS2V(B,W,N,P0,DRS) 01291 01292 P1=P0+(S2-RS)/DRS 01293 01294 IF (ABS(P1-P0)/P0.LT.EPS) THEN 01295 IF (P1.GE.0.) THEN 01296 IREP=0 01297 P=P1 01298 P0=P1 01299 ELSE 01300 IREP=-3 01301 P=0. 01302 P0=0. 01303 WRITE(*,FMT='(A34,G15.5)') & 01304 "SPLPS2V : SOLUTION NEGATIVE : P = ",P 01305 ENDIF 01306 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPS2V',1,ZHOOK_HANDLE) 01307 RETURN 01308 ENDIF 01309 ENDDO 01310 01311 IREP=-5 01312 P=P1 01313 P0=0. 01314 WRITE(*,FMT='(A48,A23,G15.5)') & 01315 "SPLPS2V : NOMBRE MAXIMAL D ITERATIONS ATTEINT : ",& 01316 "VALEUR DE P ATTEINTE : ",P1 01317 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPS2V',1,ZHOOK_HANDLE) 01318 01319 END SUBROUTINE SPLPS2V 01320 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01321 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01322 !PROCEDURES INTERMEDIAIRES POUR SPLPR 01323 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01324 !22. Comme SPLS2VI mais on ne divise pas par SIZE(W) mais par N 01325 SUBROUTINE SPLRI(W,N,RES) 01326 !SPLPR 01327 IMPLICIT NONE 01328 01329 REAL, DIMENSION(:), INTENT(IN) :: W 01330 INTEGER, INTENT(IN) :: N 01331 REAL, INTENT(OUT) :: RES 01332 01333 INTEGER :: I 01334 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01335 01336 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLRI',0,ZHOOK_HANDLE) 01337 RES=0. 01338 DO I=1,SIZE(W) 01339 RES = RES + W(I)**2 01340 ENDDO 01341 RES=RES/N 01342 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLRI',1,ZHOOK_HANDLE) 01343 01344 END SUBROUTINE SPLRI 01345 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01346 !23. RS=SQRT(S2/(SOMME((W(I)/B(I))**2)*N) 01347 SUBROUTINE SPLPR0(B,W,N,S2,P) 01348 !SPLPR 01349 IMPLICIT NONE 01350 01351 REAL, DIMENSION(:), INTENT(IN) :: B 01352 REAL, DIMENSION(:), INTENT(IN) :: W 01353 INTEGER, INTENT(IN) :: N 01354 REAL, INTENT(IN) :: S2 01355 REAL, INTENT(OUT) :: P 01356 01357 REAL :: RES 01358 INTEGER :: I 01359 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01360 01361 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPR0',0,ZHOOK_HANDLE) 01362 RES=0. 01363 DO I=1,SIZE(W) 01364 RES=RES+(W(I)/B(I))**2 01365 ENDDO 01366 P=SQRT(S2/(RES*N)) 01367 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPR0',1,ZHOOK_HANDLE) 01368 01369 END SUBROUTINE SPLPR0 01370 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01371 !24. RS=N*SOMME((W(I)/(B(I)+N*P))**2)*P**2 01372 SUBROUTINE SPLRS(B,W,N,P,RES) 01373 !SPLP, SPLPR, SP0CVQ 01374 IMPLICIT NONE 01375 01376 REAL, DIMENSION(:), INTENT(IN) :: B 01377 REAL, DIMENSION(:), INTENT(IN) :: W 01378 INTEGER, INTENT(IN) :: N 01379 REAL, INTENT(IN) :: P 01380 REAL, INTENT(OUT):: RES 01381 01382 INTEGER :: I 01383 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01384 01385 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLRS',0,ZHOOK_HANDLE) 01386 RES=0. 01387 DO I=1,SIZE(B) 01388 RES=RES+(W(I)/(B(I)+N*P))**2 01389 ENDDO 01390 RES=N*RES*P**2 01391 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLRS',1,ZHOOK_HANDLE) 01392 01393 END SUBROUTINE SPLRS 01394 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01395 !25. DRS=2*N*P*(SOMME(W(I)**2*B(I)/(B(I)+N*P)**3) 01396 SUBROUTINE SPLDRS(B,W,N,P,RES) 01397 !SPLPR 01398 IMPLICIT NONE 01399 01400 REAL, DIMENSION(:), INTENT(IN) :: B 01401 REAL, DIMENSION(:), INTENT(IN) :: W 01402 INTEGER, INTENT(IN) :: N 01403 REAL, INTENT(IN) :: P 01404 REAL, INTENT(OUT) :: RES 01405 ! 01406 INTEGER :: J 01407 REAL :: D 01408 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01409 01410 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLDRS',0,ZHOOK_HANDLE) 01411 RES=0. 01412 DO J=1,SIZE(B) 01413 D=B(J)+N*P 01414 RES=RES+W(J)**2*B(J)/(D**3) 01415 ENDDO 01416 RES=2.*N*P*RES 01417 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLDRS',1,ZHOOK_HANDLE) 01418 01419 END SUBROUTINE SPLDRS 01420 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01421 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01422 !26. SPLPR procedure: parallèle de splps2v 01423 SUBROUTINE SPLPR(N,W,B,P0,S2,P,IREP) 01424 !SPLP 01425 IMPLICIT NONE 01426 01427 INTEGER, INTENT(IN) :: N 01428 REAL, DIMENSION(:), INTENT(IN) :: W !N-M 01429 REAL, DIMENSION(:), INTENT(IN) :: B !N-M 01430 REAL, INTENT(INOUT) :: P0 01431 REAL, INTENT(INOUT) :: S2 01432 REAL, INTENT(INOUT) :: P 01433 INTEGER, INTENT(OUT) :: IREP 01434 01435 INTEGER, PARAMETER :: EPS=1.E-2 01436 REAL :: RS, DRS, P1, RINF 01437 INTEGER :: NITER 01438 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01439 01440 ! Calcul de la valeur de la fonction de Reinsch pour p=infini 01441 01442 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPR',0,ZHOOK_HANDLE) 01443 CALL SPLRI(W,N,RINF) 01444 IF (RINF.LE.S2) THEN 01445 IREP=-6 01446 P0=0. 01447 WRITE(*,FMT='(A51,G15.5)') & 01448 "SPLPR : S2 > VALEUR DE LA FONCTION POUR P INFINI = ",RINF 01449 P=(B(SIZE(B))**2)/B(1) 01450 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPR',1,ZHOOK_HANDLE) 01451 RETURN 01452 ENDIF 01453 01454 ! Recherche du point de depart de l'iteration 01455 IF (P0.EQ.0.) CALL SPLPR0(B,W,N,S2,P0) 01456 P1=P0 01457 01458 ! Iterations de Newton 01459 DO NITER=1,NSDMAX 01460 P0=P1 01461 CALL SPLRS(B,W,N,P0,RS) 01462 CALL SPLDRS(B,W,N,P0,DRS) 01463 01464 P1=P0+(S2-RS)/DRS 01465 01466 IF (ABS(P1-P0)/P0.LT.EPS) THEN 01467 IF (P1.GE.0.) THEN 01468 IREP=0 01469 P=P1 01470 P0=P1 01471 ELSE 01472 IREP=-3 01473 P=0. 01474 P0=0. 01475 WRITE(*,FMT='(A32,G15.5)') & 01476 "SPLPR : SOLUTION NEGATIVE : P = ",P 01477 ENDIF 01478 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPR',1,ZHOOK_HANDLE) 01479 RETURN 01480 ENDIF 01481 ENDDO 01482 01483 IREP=-5 01484 P=0. 01485 P0=0. 01486 WRITE(*,FMT='(A45,A23,G15.5)') & 01487 "SPLPR : NOMBRE MAXIMAL D ITERATIONS ATTEINT: ",& 01488 "VALEUR DE P ATTEINTE : ",P1 01489 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPR',1,ZHOOK_HANDLE) 01490 01491 END SUBROUTINE SPLPR 01492 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01493 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01494 !PROCEDURES INTERMEDIAIRES POUR SPLPV 01495 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01496 !27. 01497 SUBROUTINE SPLDV(B,W,N,P,RES) 01498 !SPLPV 01499 IMPLICIT NONE 01500 01501 REAL, DIMENSION(:), INTENT(IN) :: B 01502 REAL, DIMENSION(:), INTENT(IN) :: W 01503 INTEGER, INTENT(IN) :: N 01504 REAL, INTENT(IN) :: P 01505 REAL, INTENT(OUT) :: RES 01506 ! 01507 REAL :: S1, S2, SB2, SB3, A1, A2, B2 01508 INTEGER :: I 01509 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01510 01511 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLDV',0,ZHOOK_HANDLE) 01512 S1=0. 01513 S2=0. 01514 SB2=0. 01515 SB3=0. 01516 DO I=1,SIZE(B) 01517 A1=1./(B(I)+N*P) 01518 A2=A1**2 01519 B2=A2*W(I)**2 01520 S1=S1+A1 01521 S2=S2+A2 01522 SB2=SB2+B2 01523 SB3=SB3+B2*A1 01524 ENDDO 01525 RES=S2*SB2-S1*SB3 01526 RES=2.*N**2*RES/(S1**3) 01527 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLDV',1,ZHOOK_HANDLE) 01528 01529 END SUBROUTINE SPLDV 01530 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01531 !28. 01532 SUBROUTINE SPLD2V(B,W,N,P,RES) 01533 !SPLPV 01534 IMPLICIT NONE 01535 01536 REAL, DIMENSION(:), INTENT(IN) :: B 01537 REAL, DIMENSION(:), INTENT(IN) :: W 01538 INTEGER, INTENT(IN) :: N 01539 REAL, INTENT(IN) :: P 01540 REAL, INTENT(OUT) :: RES 01541 01542 REAL :: S1, S2, S3, SB2, SB3, SB4 01543 REAL :: A1, A2, B2 01544 INTEGER :: I 01545 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01546 01547 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLD2V',0,ZHOOK_HANDLE) 01548 S1=0. 01549 S2=0. 01550 S3=0. 01551 SB2=0. 01552 SB3=0. 01553 SB4=0. 01554 DO I=1,SIZE(B) 01555 A1=1./(B(I)+N*P) 01556 A2=A1**2 01557 B2=A2*W(I)**2 01558 S1=S1+A1 01559 S2=S2+A2 01560 S3=S3+A2*A1 01561 SB2=SB2+B2 01562 SB3=SB3+B2*A1 01563 SB4=SB4+B2*A2 01564 ENDDO 01565 RES=3.*S2**2*SB2-4.*S2*S1*SB3+3.*S1**2*SB4-2.*S3*S1*SB2 01566 RES=RES/(S1**4) 01567 RES=RES*2.*N**3 01568 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLD2V',1,ZHOOK_HANDLE) 01569 01570 END SUBROUTINE SPLD2V 01571 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01572 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01573 !29. 01574 SUBROUTINE SPLPV(N,B,W,P0,P,IREP) 01575 !SPLP 01576 IMPLICIT NONE 01577 01578 INTEGER, INTENT(IN) :: N 01579 REAL, DIMENSION(:), INTENT(IN):: B, W !N-M 01580 REAL, INTENT(INOUT) :: P0 01581 REAL, INTENT(OUT) :: P 01582 INTEGER, INTENT(OUT) :: IREP 01583 01584 INTEGER, PARAMETER:: EPS=1.E-100, EPSR=1.E-2 01585 REAL :: DVM, D2VM, P1 01586 INTEGER :: IFLAG, NITER 01587 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01588 01589 01590 ! Calcul de Vm'(0) . Si Vm'(0) > 0 alors popt = 0. 01591 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPV',0,ZHOOK_HANDLE) 01592 CALL SPLDV(B,W,N,0.,DVM) 01593 01594 IF (DVM.GT.0.) THEN 01595 IREP=0. 01596 P=0. 01597 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPV',1,ZHOOK_HANDLE) 01598 RETURN 01599 ENDIF 01600 01601 ! Recherche du point de depart p0 de l'iteration de Newton 01602 IFLAG=0 01603 01604 DO WHILE(IFLAG==0) 01605 01606 IF (P0.EQ.0.) THEN 01607 01608 IFLAG=1 01609 P0=B(1)/N 01610 01611 !on augmente P0 et on calcul D2VM; quand D2VM est supérieur à 0, 01612 !on passe à la suite 01613 DO NITER=1,NSDMAX 01614 CALL SPLD2V(B,W,N,P0,D2VM) 01615 IF (D2VM.GT.0.) EXIT 01616 P0=P0*10. 01617 ENDDO 01618 01619 IF (D2VM.LE.0.) THEN 01620 IREP=-1 01621 P0=0. 01622 P=0. 01623 WRITE(*,FMT='(A41)') & 01624 "SPLPV : PAS DE MINIMUM: D2VM < 0 PARTOUT " 01625 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPV',1,ZHOOK_HANDLE) 01626 RETURN 01627 ENDIF 01628 01629 ENDIF 01630 01631 P1=P0 01632 01633 DO NITER=1,NSDMAX 01634 01635 P0=P1 01636 CALL SPLDV(B,W,N,P0,DVM) 01637 CALL SPLD2V(B,W,N,P0,D2VM) 01638 01639 !Si D2Vm est très proche de 0, on sort et P=P0 01640 IF (ABS(D2VM).LT.EPS) THEN 01641 P0=0. 01642 IF (IFLAG.EQ.0) EXIT 01643 IREP=-2 01644 P=P0 01645 WRITE(*,FMT='(A28)') "SPLPV : PASSAGE PAR D2VM = 0" 01646 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPV',1,ZHOOK_HANDLE) 01647 RETURN 01648 ENDIF 01649 01650 P1=P0-DVM/D2VM 01651 01652 !si l'écart entre P0 et P1 est suffisamment petit 01653 IF (ABS(P1-P0)/P0.LT.EPSR) THEN 01654 !Si P1 est négatif, on repart dans la boucle 01655 IF (P1.LT.0.) THEN 01656 P0=0. 01657 IF (IFLAG.EQ.0) EXIT 01658 IREP=-3 01659 P=0. 01660 WRITE(*,FMT='(A33,G15.5)') & 01661 "SPLPV : SOLUTION NEGATIVE : P = ",P1 01662 ELSE 01663 !si P1 et positif et D2VM positif, on garde P1 01664 IF (D2VM.GE.0.) THEN 01665 IREP=0 01666 P=P1 01667 P0=P1 01668 ELSE 01669 !si P1 est négatif, on repart dans la boucle 01670 P0=0. 01671 IF (IFLAG.EQ.0) EXIT 01672 IREP=-4 01673 P=0. 01674 WRITE(*,FMT='(A28,G15.5)') & 01675 "SPLPV : MAXIMUM DE VM : P = ",P1 01676 ENDIF 01677 ENDIF 01678 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPV',1,ZHOOK_HANDLE) 01679 RETURN 01680 ENDIF 01681 ENDDO 01682 01683 P0=0. 01684 IF (IFLAG.NE.0) THEN 01685 IREP=-5 01686 P=P1 01687 WRITE(*,FMT='(A46,A23,G15.5)') & 01688 "SPLPV : NOMBRE MAXIMAL D ITERATIONS ATTEINT : ",& 01689 "VALEUR DE P ATTEINTE : ",P1 01690 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPV',1,ZHOOK_HANDLE) 01691 RETURN 01692 ENDIF 01693 01694 ENDDO 01695 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLPV',1,ZHOOK_HANDLE) 01696 01697 END SUBROUTINE SPLPV 01698 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01699 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01700 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01701 !30. 01702 SUBROUTINE SPLP(LISSAGE,N,B,W,P,S2,IREP) 01703 !SP0NOP 01704 IMPLICIT NONE 01705 01706 INTEGER, INTENT(IN) :: LISSAGE 01707 INTEGER, INTENT(IN) :: N 01708 REAL, DIMENSION(:), INTENT(IN) :: B 01709 REAL, DIMENSION(:), INTENT(IN) :: W 01710 REAL, INTENT(INOUT) :: P 01711 REAL, INTENT(INOUT) :: S2 01712 INTEGER, INTENT(OUT) :: IREP 01713 01714 REAL :: P0 01715 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01716 01717 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLP',0,ZHOOK_HANDLE) 01718 IF (LISSAGE.EQ.1) THEN 01719 ! Recherche de la valeur de P pour S2 connu, par la methode de Wahba 01720 P0=0. 01721 CALL SPLPS2V(B,W,N,P0,S2,P,IREP) 01722 IF (IREP.NE.0) THEN 01723 IREP=-2 01724 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLP',1,ZHOOK_HANDLE) 01725 RETURN 01726 ENDIF 01727 ELSEIF (LISSAGE.EQ.10) THEN 01728 ! Recherche de la valeur de P pour S2 connu, par la methode de Reinsch 01729 P0=0. 01730 CALL SPLPR(N,W,B,P0,S2,P,IREP) 01731 IF (IREP.NE.0) THEN 01732 IREP=-2 01733 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLP',1,ZHOOK_HANDLE) 01734 RETURN 01735 ENDIF 01736 ELSEIF (LISSAGE.EQ.2) THEN 01737 ! Recherche de la valeur de P par la methode de G.C.V 01738 P0=0. 01739 CALL SPLPV(N,B,W,P0,P,IREP) 01740 IF (IREP.NE.0) IREP=-3 01741 ! Calcul de la variance d'erreur de mesure correspondante 01742 CALL SPLS2V(B,W,N,P,S2) 01743 ELSEIF (LISSAGE.EQ.3) THEN 01744 ! Calcul de l'ajustement correspondant a une valeur de P donnee 01745 CALL SPLRS(B,W,N,P,S2) 01746 ENDIF 01747 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLP',1,ZHOOK_HANDLE) 01748 01749 END SUBROUTINE SPLP 01750 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01751 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01752 ! GROUPE DES PROCEDURES T et R 01753 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01754 !31. TT procedure 01755 SUBROUTINE SPLTT(X,DS,IW,TM,TT,IREP) 01756 !SPLR 01757 IMPLICIT NONE 01758 01759 REAL,DIMENSION(:,:),INTENT(IN) :: X !ND,N 01760 REAL,DIMENSION(:),INTENT(IN) :: DS !N 01761 INTEGER,DIMENSION(:,:),INTENT(IN) :: IW !ND,M 01762 REAL,DIMENSION(:,:),INTENT(OUT) :: TM !N,M 01763 REAL,DIMENSION(:,:),INTENT(OUT) :: TT !M,M 01764 INTEGER, INTENT(OUT) :: IREP 01765 01766 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01767 01768 ! Calcul de la matrice T des monomes 01769 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLTT',0,ZHOOK_HANDLE) 01770 CALL SPLT(X,IW,TM) 01771 ! Calcul de Tm=(Ds-1)*T 01772 CALL MXIDML(DS,TM,TM) 01773 ! Calcul de Tm'*Tm 01774 CALL MXMSPL(TRANSPOSE(TM),TM,TT) 01775 ! Calcul de (Tm'*Tm)-1 01776 CALL SMXINV(TT,IREP) 01777 ! Test du compte-rendu de smxinv 01778 IF (IREP.NE.0) WRITE(*,FMT='(A27)') "SPLTT: MATRICE TT NON INVERSIBLE" 01779 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLTT',1,ZHOOK_HANDLE) 01780 01781 END SUBROUTINE SPLTT 01782 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01783 !32. r procedure 01784 SUBROUTINE SPLR_1(NORD,X,DS,IW,TTT,R,IREP) 01785 !SP0NOP, SP0CVQ 01786 IMPLICIT NONE 01787 01788 INTEGER,INTENT(IN) :: NORD 01789 REAL,DIMENSION(:,:), INTENT(IN) :: X !ND,N 01790 REAL,DIMENSION(:),INTENT(IN) :: DS !N 01791 INTEGER,DIMENSION(:),INTENT(OUT) :: IW !ND 01792 REAL,DIMENSION(:),INTENT(OUT) :: TTT !N 01793 REAL,DIMENSION(:,:),INTENT(OUT) :: R !N,N 01794 INTEGER, INTENT(OUT)::IREP 01795 01796 INTEGER :: N, I, J 01797 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01798 01799 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLR_1',0,ZHOOK_HANDLE) 01800 N=SIZE(X,2) 01801 01802 ! R = I 01803 DO I=1,N 01804 DO J=1,N 01805 IF (I.EQ.J) THEN 01806 R(I,I)=1. 01807 ELSE 01808 R(I,J)=0. 01809 ENDIF 01810 ENDDO 01811 ENDDO 01812 01813 IREP=0 01814 01815 ! Calcul de (Ds-1)*R 01816 CALL MXIDML(DS,R,R) 01817 01818 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLR_1',1,ZHOOK_HANDLE) 01819 01820 END SUBROUTINE SPLR_1 01821 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01822 !32. 01823 SUBROUTINE SPLR_2(NORD,X,DS,IW,TTT,R,IREP) 01824 01825 IMPLICIT NONE 01826 01827 INTEGER,INTENT(IN) :: NORD 01828 REAL,DIMENSION(:,:), INTENT(IN) :: X !ND,N 01829 REAL,DIMENSION(:),INTENT(IN) :: DS !N 01830 INTEGER,DIMENSION(:,:),INTENT(OUT) :: IW !ND,M 01831 REAL,DIMENSION(:,:),INTENT(OUT) :: TTT !M,N 01832 REAL,DIMENSION(:,:),INTENT(OUT) :: R !N,N-M 01833 INTEGER, INTENT(OUT)::IREP 01834 01835 REAL,DIMENSION(SIZE(X,2),SIZE(X,2)) :: C !N,N 01836 REAL,DIMENSION(SIZE(X,2),SIZE(IW,2)) :: TM !N,M 01837 REAL,DIMENSION(SIZE(IW,2),SIZE(IW,2)) :: TT !M,M 01838 REAL, DIMENSION(SIZE(X,2),SIZE(X,2)) :: R1 01839 REAL,DIMENSION(SIZE(X,2)) :: VP !N 01840 REAL,DIMENSION(SIZE(X,2)) :: WORK !N 01841 INTEGER :: M, N, J 01842 INTEGER, DIMENSION(1) :: IMIN 01843 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01844 01845 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLR_2',0,ZHOOK_HANDLE) 01846 N=SIZE(X,2) 01847 M=SIZE(IW,2) 01848 01849 ! Recherche des monomes de degre inferieur ou egal a nord-1 01850 CALL SPLIE(NORD-1,IW) 01851 !---------------------------------------------------------- 01852 ! Calcul de la matrice (Tm'*Tm)-1 01853 CALL SPLTT(X,DS,IW,TM,TT,IREP) 01854 ! Calcul de (T'm*Tm)-1 Tm' 01855 CALL MXMSPL(TT,TRANSPOSE(TM),TTT) 01856 !R1(:,1:M)=TRANSPOSE(TTT) 01857 !------------------------------------------------------------ 01858 ! Calcul de Tm ((Tm'*Tm)-1) Tm' 01859 CALL MTXAXM(TRANSPOSE(TM),TT,C) 01860 DO J=1,N 01861 C(J,J)=C(J,J)-1. 01862 ENDDO 01863 C(:,:)=-C(:,:) 01864 !----------------------------------------------------------- 01865 ! Calcul des vecteurs propres de I - Tm c ((Tm'*Tm)-1) c Tm' 01866 CALL EISRS1(C,VP,R1,WORK,IREP) 01867 01868 ! Test du signe des valeurs propres de C 01869 IMIN=MINLOC(VP(M+1:N)) 01870 IF(VP(M+IMIN(1)).LE.0) THEN 01871 IREP=-1 01872 WRITE(*,FMT='(A31)') "SPLR: VALEUR PROPRE < 0 DE PROJ" 01873 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLR_2',1,ZHOOK_HANDLE) 01874 RETURN 01875 ENDIF 01876 ! Compte-rendu OK 01877 IREP=0 01878 ! Calcul de (Ds-1)*R 01879 R=R1(:,M+1:N) 01880 CALL MXIDML(DS,R,R) 01881 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLR_2',1,ZHOOK_HANDLE) 01882 01883 END SUBROUTINE SPLR_2 01884 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01885 ! GROUPE DES PROCEDURES U 01886 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01887 !33. U procedure 01888 SUBROUTINE SPLU(NORD,X,G,R,RK,U,DB) 01889 !SP0NOP, SP0CVQ 01890 IMPLICIT NONE 01891 01892 INTEGER, INTENT(IN) :: NORD 01893 REAL,DIMENSION(:,:), INTENT(IN) :: X !ND,N 01894 REAL,DIMENSION(:,:), INTENT(IN) :: G ! ND,ND 01895 REAL,DIMENSION(:,:), INTENT(IN) :: R !N,N-M 01896 REAL,DIMENSION(:,:), INTENT(OUT) :: RK !N,N 01897 REAL,DIMENSION(:,:), INTENT(OUT) :: U !N-M,N-M 01898 REAL,DIMENSION(:), INTENT(OUT) :: DB !N-M 01899 01900 REAL,DIMENSION(SIZE(R,2),SIZE(R,2)) :: RKR !N-M,N-M 01901 REAL,DIMENSION(SIZE(R,2)) :: WORK !N-M 01902 01903 INTEGER :: M, N, IREP 01904 INTEGER, DIMENSION(1) :: IMIN 01905 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01906 01907 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLU',0,ZHOOK_HANDLE) 01908 N=SIZE(X,2) 01909 M=SIZE(X,2)-SIZE(R,2) 01910 01911 ! Calcul de la matrice K des noyaux 01912 CALL SPLK(NORD,X,X,G,RK) 01913 ! Calcul de R'*K*R 01914 CALL MTXAXM(R,RK,RKR) 01915 ! Calcul de la matrice U des vecteurs propres de R'*K*R 01916 CALL EISRS1(RKR,DB,U,WORK,IREP) 01917 ! Test du signe des valeurs propres de R'KR 01918 IMIN=MINLOC(DB) 01919 IF (DB(IMIN(1)).LE.0) & 01920 WRITE(*,FMT='(A31)') "SPLU : VALEUR PROPRE < 0 DE RKR" 01921 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLU',1,ZHOOK_HANDLE) 01922 01923 END SUBROUTINE SPLU 01924 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01925 ! GROUPE DES PROCEDURES W 01926 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01927 !34. w procedure 01928 SUBROUTINE SPLW(Z,R,U,W,RZ) 01929 !SP0NOP, SP0CVQ 01930 IMPLICIT NONE 01931 01932 REAL,DIMENSION(:), INTENT(IN) :: Z !N 01933 REAL, DIMENSION(:,:),INTENT(IN) :: R !N,N-M 01934 REAL,DIMENSION(:,:),INTENT(IN) :: U !N-M,N-M 01935 REAL, DIMENSION(:),INTENT(OUT) :: W !N-M 01936 REAL,DIMENSION(:),INTENT(OUT) :: RZ !N-M 01937 01938 INTEGER :: N, M 01939 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01940 01941 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLW',0,ZHOOK_HANDLE) 01942 N=SIZE(Z) 01943 M=SIZE(Z)-SIZE(R,2) 01944 01945 ! Calcul de R'*Z 01946 CALL MXMSPL(TRANSPOSE(R),Z,RZ) 01947 ! Calcul de U'R'Z 01948 CALL MXMSPL(TRANSPOSE(U),RZ,W) 01949 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLW',1,ZHOOK_HANDLE) 01950 01951 END SUBROUTINE SPLW 01952 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01953 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01954 ! GROUPE DES PROCEDURES Q 01955 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01956 !35. splvpQ procedure 01957 SUBROUTINE SPLVPQ(P,R,RK,U,DB,Q) 01958 !SP0CVQ 01959 IMPLICIT NONE 01960 01961 REAL, INTENT(IN) :: P 01962 REAL, DIMENSION(:,:), INTENT(IN) :: R !N,N-M 01963 REAL, DIMENSION(:,:), INTENT(IN) :: RK ! N,N 01964 REAL, DIMENSION(:,:), INTENT(IN) :: U !N-M,N-M 01965 REAL, DIMENSION(:), INTENT(IN) :: DB !N-M 01966 REAL, DIMENSION(:,:), INTENT(OUT) :: Q ! N,N 01967 01968 REAL, DIMENSION(SIZE(R,1),SIZE(R,1)):: Q1 !N,N 01969 REAL, DIMENSION(SIZE(U,1),SIZE(U,2))::QW1,QW2 !N-M,N-M 01970 INTEGER :: N, J 01971 REAL(KIND=JPRB) :: ZHOOK_HANDLE 01972 01973 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLVPQ',0,ZHOOK_HANDLE) 01974 N=SIZE(R,1) 01975 01976 ! Calcul de (DB+n*p*I)-1 01977 QW1(:,:)=0. 01978 DO J=1,SIZE(DB) 01979 QW1(J,J)=1./(DB(J)+N*P) 01980 ENDDO 01981 01982 ! Calcul de U((DB+n*p*I)-1)U' 01983 CALL MTXAXM(TRANSPOSE(U),QW1,QW2) 01984 ! Calcul de RU((DB+n*p*I)-1))UR' 01985 CALL MTXAXM(TRANSPOSE(R),QW2,Q) 01986 ! Calcul de QK 01987 CALL MXMSPL(Q,RK,Q1) 01988 01989 ! Calcul de I - QK 01990 Q1(:,:)=-Q1(:,:) 01991 DO J=1,N 01992 Q1(J,J)=Q1(J,J)+1 01993 ENDDO 01994 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLVPQ',1,ZHOOK_HANDLE) 01995 01996 END SUBROUTINE SPLVPQ 01997 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 01998 ! GROUPE DES PROCEDURES C 01999 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02000 !36. 02001 SUBROUTINE SPLC(Q,Z,RK,DS,TTT,RKC,C) 02002 !SP0NOP, SP0CVQ 02003 IMPLICIT NONE 02004 02005 REAL, DIMENSION(:,:), INTENT(IN) :: Q !Q(N,N) 02006 REAL, DIMENSION(:), INTENT(IN) :: Z !Z(N) 02007 REAL, DIMENSION(:,:), INTENT(IN) :: RK !RK(N,n) 02008 REAL, DIMENSION(:), INTENT(IN) :: DS !DS(N) 02009 REAL, DIMENSION(:,:), INTENT(IN) :: TTT !TTT(M,N) 02010 REAL, DIMENSION(:), INTENT(OUT) :: RKC !RKC(N) 02011 REAL, DIMENSION(:), INTENT(OUT) :: C !C(N+M) 02012 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02013 02014 INTEGER :: N, M 02015 02016 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLC',0,ZHOOK_HANDLE) 02017 N=SIZE(Z) 02018 M=SIZE(TTT,1) 02019 02020 ! Calcul de c 02021 ! ----------- 02022 ! Calcul de c = Qz = R((R'KR+n*p*I)-1)R'z 02023 CALL MXMSPL(Q,Z,C(1:N)) 02024 02025 ! Calcul de d 02026 ! ----------- 02027 IF (M.NE.0) THEN 02028 ! Calcul de K*c 02029 CALL MXMSPL(RK,C(1:N),RKC) 02030 ! Calcul de z - K*c 02031 RKC(:)=Z(:)-RKC(:) 02032 ! Calcul de (Ds-1) * (z-K*c) 02033 RKC(:)=RKC(:)/DS(:) 02034 ! Calcul de d = (Tm'*Tm)-1 * Tm' * (Ds-1) * (z-K*c) 02035 CALL MXMSPL(TTT,RKC,C(N+1:N+M)) 02036 ENDIF 02037 02038 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLC',1,ZHOOK_HANDLE) 02039 02040 END SUBROUTINE SPLC 02041 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02042 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02043 ! PROCEDURES GLOBALES 02044 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02045 !37. Calcul de M en fonction de NORD et ND (qu'est-ce que c'est que ce ND?) 02046 SUBROUTINE SPLM(ND,NORD,M) 02047 !SP0NOP 02048 IMPLICIT NONE 02049 02050 INTEGER, INTENT(IN) :: ND 02051 INTEGER, INTENT(IN) :: NORD 02052 INTEGER, INTENT(OUT):: M 02053 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02054 02055 02056 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLM',0,ZHOOK_HANDLE) 02057 IF (ND==1) THEN 02058 M=NORD 02059 ELSEIF (ND==2) THEN 02060 M=NORD*(NORD+1)/2. 02061 ELSEIF (ND==3) THEN 02062 M=NORD*(NORD+1)*(NORD+2)/6. 02063 ELSEIF (ND==4) THEN 02064 M=NORD*(NORD+1)*(NORD+2)*(NORD+3)/24. 02065 ENDIF 02066 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLM',1,ZHOOK_HANDLE) 02067 02068 END SUBROUTINE SPLM 02069 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02070 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02071 !38. 02072 SUBROUTINE SP0NOP(X,G,Z,DS,LISSAGE,LORDRE,IOPT,NORDOPT,M,S2,P,& 02073 IW,TTT,R,RK,U,DB,W,RZ,IREP) 02074 !SP0CVQ 02075 IMPLICIT NONE 02076 02077 REAL, DIMENSION(:,:), INTENT(IN) :: X !ND,N 02078 REAL, DIMENSION(:,:), INTENT(IN) :: G ! ND,ND 02079 REAL, DIMENSION(:), INTENT(IN) :: Z ! N 02080 REAL, DIMENSION(:), INTENT(IN) :: DS ! N 02081 INTEGER, INTENT(IN) :: LISSAGE 02082 INTEGER, INTENT(IN) :: LORDRE 02083 INTEGER, INTENT(IN) :: IOPT 02084 INTEGER, INTENT(INOUT) :: NORDOPT 02085 INTEGER, INTENT(INOUT) :: M 02086 REAL, INTENT(INOUT) :: S2 02087 REAL, INTENT(INOUT) :: P 02088 INTEGER, DIMENSION(:,:), INTENT(OUT):: IW !ND,M 02089 REAL, DIMENSION(:,:), INTENT(OUT) :: TTT ! M,N 02090 REAL, DIMENSION(:,:), INTENT(OUT) :: R ! N,N-M 02091 REAL,DIMENSION(:,:), INTENT(OUT) :: RK !N,N 02092 REAL,DIMENSION(:,:), INTENT(OUT) :: U !N-M,N-M 02093 REAL,DIMENSION(:), INTENT(OUT) :: DB !N-M 02094 REAL, DIMENSION(:),INTENT(OUT) :: W !N-M 02095 REAL,DIMENSION(:),INTENT(OUT) :: RZ !N-M 02096 INTEGER, INTENT (OUT) :: IREP 02097 02098 02099 REAL, DIMENSION(NORDMAX) :: VM, S2SAVE, PSAVE, MSAVE 02100 INTEGER :: N, ND 02101 INTEGER :: NORDMI, NORDMA, NORDM, NORD 02102 REAL :: VMMIN 02103 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02104 02105 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SP0NOP',0,ZHOOK_HANDLE) 02106 N=SIZE(X,2) 02107 ND=SIZE(X,1) 02108 02109 IF (LORDRE.EQ.0) THEN 02110 ! 02111 ! Recherche de p seulement avec nord fixe 02112 ! --------------------------------------- 02113 IF (IOPT.EQ.0) THEN 02114 ! Calcul de IW, TTT, R 02115 CALL SPLR(NORDOPT,X,DS,IW,TTT,R,IREP) 02116 ! Calcul de RK, U, DB 02117 CALL SPLU(NORDOPT,X,G,R,RK,U,DB) 02118 ENDIF 02119 !calcul de W, RZ 02120 CALL SPLW(Z,R,U,W,RZ) 02121 !calcul de P 02122 CALL SPLP(LISSAGE,N,DB,W,P,S2,IREP) 02123 02124 ELSE 02125 02126 ! Recherche de nordopt et popt 02127 ! ---------------------------- 02128 02129 NORDMI=ND/2+1 02130 02131 ! Calcul de l'ordre maximal pour le nombre de points n 02132 CALL SPLM(ND,NORD,M) 02133 DO WHILE(N.GT.M) 02134 NORD=NORD+1 02135 CALL SPLM(ND,NORD,M) 02136 ENDDO 02137 NORDMA=MIN(NORD,NORDMAX) 02138 NORDM=NORDMA 02139 02140 ! Calcul de Vm pour les differents ordre 02141 DO NORD=NORDMI,NORDMA 02142 !calcul de M 02143 CALL SPLM(ND,NORD,M) 02144 CALL SPLR(NORD,X,DS,IW,TTT,R,IREP) 02145 CALL SPLU(NORD,X,G,R,RK,U,DB) 02146 CALL SPLW(Z,R,U,W,RZ) 02147 CALL SPLP(LISSAGE,N,DB,W,P,S2,IREP) 02148 MSAVE(NORD)=M 02149 S2SAVE(NORD)=S2 02150 PSAVE(NORD)=P 02151 IF (IREP.EQ.0) THEN 02152 !calcul de VM(NORD) 02153 CALL SPLV(DB,W,N,P,VM(NORD)) 02154 WRITE(*,FMT='(A15,I5,A5,G15.5,A6,G15.5)')& 02155 "CNORD : NORD = ",NORD," P = ",P," VM = ",VM(NORD) 02156 ELSE 02157 IF (NORD.EQ.NORDMI) THEN 02158 NORDOPT=NORD 02159 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SP0NOP',1,ZHOOK_HANDLE) 02160 RETURN 02161 ELSE 02162 NORDM=NORD-1 02163 IREP=0 02164 EXIT 02165 ENDIF 02166 ENDIF 02167 ENDDO 02168 02169 ! Recherche de l'ordre optimal 02170 VMMIN=1.E200 02171 DO NORD=NORDMI,NORDM 02172 IF(VM(NORD).LT.VMMIN) THEN 02173 VMMIN=VM(NORD) 02174 NORDOPT=NORD 02175 M=MSAVE(NORD) 02176 S2=S2SAVE(NORD) 02177 P=PSAVE(NORD) 02178 ENDIF 02179 ENDDO 02180 WRITE(*,FMT='(A24,I5)') "CNORD : ORDRE OPTIMAL = ",NORDOPT 02181 ENDIF 02182 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SP0NOP',1,ZHOOK_HANDLE) 02183 02184 END SUBROUTINE SP0NOP 02185 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02186 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02187 !39. 02188 SUBROUTINE SP0CVQ(NORD,M,X,G,Z,DS,S2,P,LISSAGE,LORDRE,IOPT,C,IREP) 02189 !SPLB2C 02190 IMPLICIT NONE 02191 02192 INTEGER, INTENT(INOUT) :: NORD 02193 INTEGER, INTENT(INOUT) :: M 02194 REAL,DIMENSION(:,:), INTENT(IN) :: X !ND,N 02195 REAL, DIMENSION(:,:), INTENT(IN) :: G !ND,ND 02196 REAL,DIMENSION(:), INTENT(IN) :: Z ! N 02197 REAL, DIMENSION(:), INTENT(IN) :: DS !N 02198 REAL, INTENT(INOUT) :: S2 02199 REAL, INTENT(INOUT) :: P 02200 INTEGER, INTENT(IN) :: LISSAGE 02201 INTEGER, INTENT(IN) :: LORDRE 02202 INTEGER, INTENT(IN) :: IOPT 02203 REAL, DIMENSION(:), INTENT(OUT) :: C !N+M 02204 INTEGER, INTENT(OUT) :: IREP 02205 02206 INTEGER, DIMENSION(SIZE(X,1),SIZE(C)-SIZE(X,2)) :: IW 02207 REAL, DIMENSION(SIZE(C)-SIZE(X,2),SIZE(X,2)) :: TTT ! M,N 02208 REAL, DIMENSION(SIZE(X,2),2*SIZE(X,2)-SIZE(C)) :: R ! N,N-M 02209 REAL,DIMENSION(SIZE(X,2),SIZE(X,2)) :: RK !N,N 02210 REAL,DIMENSION(2*SIZE(X,2)-SIZE(C),2*SIZE(X,2)-SIZE(C)) :: U !N-M,N-M 02211 REAL,DIMENSION(2*SIZE(X,2)-SIZE(C)):: DB !N-M 02212 REAL, DIMENSION(2*SIZE(X,2)-SIZE(C)) :: W !N-M 02213 REAL,DIMENSION(SIZE(R,2)) :: RZ !N-M 02214 REAL, DIMENSION(SIZE(X,2),SIZE(X,2)) :: Q ! N,N 02215 REAL :: RN 02216 REAL, DIMENSION(SIZE(X,2)) :: RKC !RKC(N) 02217 02218 02219 INTEGER :: N 02220 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02221 02222 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SP0CVQ',0,ZHOOK_HANDLE) 02223 N=SIZE(X,2) 02224 02225 ! Test des differents parametres 02226 IF (NORD.LT.1) THEN 02227 WRITE(*,FMT='(A17)') "SP0CVQ : NORD < 1" 02228 ELSEIF (M.LT.0) THEN 02229 WRITE(*,'(A14)') "SP0CVQ : M < 0" 02230 ELSEIF (N.LT.M) THEN 02231 WRITE(*,'(A14)') "SP0CVQ : N < M" 02232 ELSEIF (S2.LT.0.) THEN 02233 WRITE(*,'(A15)') "SP0CVQ : S2 < 0" 02234 ELSEIF (P.LT.0.) THEN 02235 WRITE(*,'(A14)') "SP0CVQ : P < 0" 02236 ELSEIF ((LISSAGE.LT.0.OR.LISSAGE.GT.3).AND.LISSAGE.NE.10) THEN 02237 WRITE(*,'(A27)') "SP0CVQ : LISSAGE < 0 OU > 3" 02238 ELSEIF (LORDRE.LT.0.OR.LORDRE.GT.1) THEN 02239 WRITE(*,'(A26)') "SP0CVQ : LORDRE < 0 OU > 1" 02240 ELSEIF (IOPT.LT.0.OR.IOPT.GT.1) THEN 02241 WRITE(*,'(A24)') "SP0CVQ : IOPT < 0 OU > 1" 02242 ELSE 02243 02244 ! Determination de nord opt ou de p pot 02245 ! ----------------------------------------- 02246 IF (LORDRE.EQ.1 .OR. & 02247 LISSAGE.EQ.1 .OR. LISSAGE.EQ.2 .OR. LISSAGE.EQ.10) THEN 02248 !calcul de P en passant par SPLR, SPLU, SPLW, si LORDRE = 0 02249 !calcul de NORD, M, S2, P si LORDRE=1 02250 CALL SP0NOP(X,G,Z,DS,LISSAGE,LORDRE,IOPT,NORD,M,S2,P,& 02251 IW,TTT,R,RK,U,DB,W,RZ,IREP) 02252 ELSEIF (LISSAGE.EQ.0) THEN 02253 P=0. 02254 ENDIF 02255 !--------------------------------------------------------------------- 02256 !--------------------------------------------------------------------- 02257 02258 IF (LORDRE.EQ.1 .OR. LISSAGE.EQ.0 .OR. LISSAGE.EQ.3) THEN 02259 IF (IOPT.EQ.0) THEN 02260 ! Calcul de IW, TTT, R 02261 CALL SPLR(NORD,X,DS,IW,TTT,R,IREP) 02262 IF (IREP.NE.0 .AND. LHOOK) CALL DR_HOOK('MODE_SPLINES:SP0CVQ',1,ZHOOK_HANDLE) 02263 IF (IREP.NE.0) RETURN 02264 ! Calcul de RK, U, DB 02265 CALL SPLU(NORD,X,G,R,RK,U,DB) 02266 ENDIF 02267 !--------------------------------------------------------------------- 02268 ! Calcul de W,RZ 02269 CALL SPLW(Z,R,U,W,RZ) 02270 ! Calcul de S2 pour lissage=3 02271 IF (LISSAGE.EQ.3) CALL SPLRS(DB,W,N,P,S2) 02272 ENDIF 02273 02274 !--------------------------------------------------------------------- 02275 !--------------------------------------------------------------------- 02276 ! Calcul de Q 02277 CALL SPLVPQ(P,R,RK,U,DB,Q) 02278 ! Calcul de RN 02279 CALL SPLBVM(DB,W,N,P,RN) 02280 ! Calcul des coefficients C et RKC 02281 CALL SPLC(Q,Z,RK,DS,TTT,RKC,C) 02282 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SP0CVQ',1,ZHOOK_HANDLE) 02283 RETURN 02284 02285 ENDIF 02286 02287 IREP=-1 02288 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SP0CVQ',1,ZHOOK_HANDLE) 02289 02290 END SUBROUTINE SP0CVQ 02291 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02292 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02293 !40. AI(ND,ND,NSDMAX): calcul de AI en fonction de XD, NSD, INTER; 02294 !remarque: on dirait que la deuxième dimension de AI est plutôt 2 que ND 02295 SUBROUTINE SPLBSD(NSD,INTER,XD,AI) 02296 !SPLB2C 02297 IMPLICIT NONE 02298 02299 INTEGER, DIMENSION(:), INTENT(IN) :: NSD 02300 INTEGER, INTENT(IN) :: INTER 02301 REAL, DIMENSION(:,:), INTENT(IN) :: XD 02302 REAL, DIMENSION(:,:,:), INTENT(OUT) :: AI 02303 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02304 02305 REAL, DIMENSION(SIZE(NSD)) :: DXI, DXR 02306 INTEGER :: J, I 02307 02308 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBSD',0,ZHOOK_HANDLE) 02309 DO J=1,SIZE(NSD) 02310 DXI(J)=(XD(2,J)-XD(1,J))/NSD(J) 02311 DXR(J)=DXI(J)/(2*INTER-2) 02312 ENDDO 02313 02314 DO J=1,SIZE(NSD) 02315 AI(J,1,1)=XD(1,J) 02316 DO I=2,NSD(J) 02317 AI(J,1,I)=XD(1,J)+(I-1)*DXI(J)-DXR(J) 02318 ENDDO 02319 AI(J,2,NSD(J))=XD(2,J) 02320 DO I=1,NSD(J)-1 02321 AI(J,2,I)=XD(2,J)-(NSD(J)-I)*DXI(J)+DXR(J) 02322 ENDDO 02323 ENDDO 02324 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBSD',1,ZHOOK_HANDLE) 02325 02326 END SUBROUTINE SPLBSD 02327 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02328 !41. récupère les X et Z tels que les valeurs soient toutes comprises entre 02329 !A et B, dans XS et ZS. NS est la taille finale des tableaux. 02330 SUBROUTINE SPLBSEL(X,Z,AI,BI,NS,XS,ZS,IREP) 02331 !SPLB2C 02332 IMPLICIT NONE 02333 02334 REAL, DIMENSION(:,:), INTENT(IN) :: X !ND,N 02335 REAL, DIMENSION(:), INTENT(IN) :: Z !N 02336 REAL, DIMENSION(:), INTENT(IN) :: AI !ND 02337 REAL, DIMENSION(:), INTENT(IN) :: BI !ND 02338 INTEGER, INTENT(OUT) :: NS 02339 REAL, DIMENSION(:,:), INTENT(OUT) :: XS !ND,NMAX 02340 REAL, DIMENSION(:), INTENT(OUT) :: ZS !NMAX 02341 INTEGER, INTENT(OUT) :: IREP 02342 02343 INTEGER :: I, J, I0 02344 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02345 02346 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBSEL',0,ZHOOK_HANDLE) 02347 IREP=0 02348 NS=0 02349 DO J=1,SIZE(X,2) 02350 I0=0 02351 DO I=1,SIZE(X,1) 02352 !il faut que tous les X(I,J) soient compris entre AI et BI 02353 IF (X(I,J).LT.AI(I).OR.X(I,J).GT.BI(I)) THEN 02354 I0=1 02355 EXIT 02356 ENDIF 02357 ENDDO 02358 IF (I0.NE.1) THEN 02359 NS=NS+1 02360 IF (NS.GT.SIZE(ZS)) THEN 02361 IREP=-1 02362 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBSEL',1,ZHOOK_HANDLE) 02363 RETURN 02364 ENDIF 02365 DO I=1,SIZE(X,1) 02366 XS(I,NS)=X(I,J) 02367 ENDDO 02368 ZS(NS)=Z(J) 02369 ENDIF 02370 ENDDO 02371 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLBSEL',1,ZHOOK_HANDLE) 02372 02373 END SUBROUTINE SPLBSEL 02374 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 02375 !42. appelée par INTERPOL_SPLINES 02376 SUBROUTINE SPLB2C(NORD,M,X,G,Z,S2,P,LISSAGE,IOPT,NSDI,NSDJ,& 02377 INTER,XD,C,IREP) 02378 !INTERPOL_SPLINES 02379 IMPLICIT NONE 02380 02381 INTEGER, INTENT(INOUT) :: NORD 02382 INTEGER, INTENT(INOUT) :: M 02383 REAL, DIMENSION(:,:), INTENT(IN) :: X !ND,N 02384 REAL, DIMENSION(:,:), INTENT(IN) :: G ! ND,ND 02385 REAL, DIMENSION(:), INTENT(IN) :: Z !N 02386 REAL, INTENT(INOUT) :: S2 02387 REAL, INTENT(INOUT) :: P 02388 INTEGER, INTENT(IN) :: LISSAGE 02389 INTEGER, INTENT(IN) :: IOPT 02390 INTEGER, INTENT(IN) :: NSDI 02391 INTEGER, INTENT(IN) :: NSDJ 02392 INTEGER, INTENT(IN) :: INTER 02393 REAL, DIMENSION(:,:), INTENT(IN) :: XD !2,ND 02394 REAL, DIMENSION(:,:,:), INTENT(OUT) :: C !nmax+mmax,nsdi,nsdj 02395 INTEGER, INTENT(OUT) :: IREP 02396 ! 02397 INTEGER :: NN, ND 02398 INTEGER :: LORDRE, I, J 02399 INTEGER, DIMENSION(2) :: NSD 02400 REAL, DIMENSION(NMAX) :: DS 02401 REAL, DIMENSION(2,2) :: CI 02402 02403 REAL(KIND=JPRB) :: ZHOOK_HANDLE 02404 02405 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLB2C',0,ZHOOK_HANDLE) 02406 ND=SIZE(X,1) 02407 IF (.NOT.ALLOCATED(AI)) ALLOCATE(AI(ND,2,NSDMAX)) 02408 IF (.NOT.ALLOCATED(NS)) ALLOCATE(NS(NSDMAX,NSDMAX)) 02409 IF (.NOT.ALLOCATED(ZS)) ALLOCATE(ZS(NMAX,NSDMAX,NSDMAX)) 02410 IF (.NOT.ALLOCATED(XS)) ALLOCATE(XS(ND,NMAX,NSDMAX,NSDMAX)) 02411 02412 ! Initialisation de Ds 02413 DS(:)=1. 02414 02415 ! Generation des sous-domaines: matrice AI de MODD_SPLINES 02416 NSD(1)=NSDI 02417 NSD(2)=NSDJ 02418 CALL SPLBSD(NSD,INTER,XD,AI) 02419 02420 LORDRE=0 02421 02422 ! Calcul des splines par sous-domaines 02423 DO I=1,NSDI 02424 DO J=1,NSDJ 02425 CI(1,1)=AI(1,1,I) 02426 CI(1,2)=AI(1,2,I) 02427 CI(2,1)=AI(2,1,J) 02428 CI(2,2)=AI(2,2,J) 02429 !selection des points 02430 CALL SPLBSEL(X,Z,CI(:,1),CI(:,2),NS(I,J),XS(:,:,I,J),ZS(:,I,J),IREP) 02431 IF (IREP.NE.0) THEN 02432 WRITE(*,FMT='(A35)') "SPLB2C: NOMBRE DE POINTS TROP GRAND" 02433 IREP=-4 02434 CALL ABOR1_SFX('SPLINES: SPLB2C: NOMBRE DE POINTS TROP GRAND') 02435 ELSEIF (NS(I,J).LE.M) THEN 02436 WRITE(*,FMT='(A42)') & 02437 "SPB2C : NB DE PTS <= M : COEF NON CALCULES" 02438 ENDIF 02439 !calcul de la spline 02440 IF (NS(I,J).GT.M) THEN 02441 NN=NS(I,J) 02442 !X(ND,NN) ZS(NN) C(NN+M) IW(ND,M) W(7*NN*NN) 02443 CALL SP0CVQ(NORD,M,XS(:,1:NN,I,J),G,ZS(1:NN,I,J),DS(1:NN),S2,P,& 02444 LISSAGE,LORDRE,IOPT,C(1:NN+M,I,J),IREP) 02445 ENDIF 02446 ENDDO 02447 ENDDO 02448 IF (LHOOK) CALL DR_HOOK('MODE_SPLINES:SPLB2C',1,ZHOOK_HANDLE) 02449 02450 END SUBROUTINE SPLB2C 02451 02452 02453 END MODULE MODE_SPLINES
1.8.0
