berkho.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\artemis\berkho.f
00002 !
00259                      SUBROUTINE BERKHO
00260 !                    *****************
00261 !
00262      &(LF)
00263 !
00264 !***********************************************************************
00265 ! ARTEMIS   V6P1                                   31/05/2011
00266 !***********************************************************************
00267 !
00268 !
00269 !
00270 !
00271 !
00272 !
00273 !
00274 !
00275 !
00276 !
00277 !
00278 !
00279 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00280 !| LF             |-->| INDICATOR OF FIRST RESOLUTION FOR DISSIPATION LOOP
00281 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00282 !
00283       USE BIEF
00284       USE DECLARATIONS_TELEMAC
00285       USE DECLARATIONS_ARTEMIS
00286       USE INTERFACE_ARTEMIS, EX_BERKHO => BERKHO
00287 !
00288       IMPLICIT NONE
00289       INTEGER LNG,LU
00290       COMMON/INFO/LNG,LU
00291 !
00292       INTEGER, INTENT(IN) :: LF
00293 !
00294       INTEGER I,ITERMU,IG
00295 
00296       DOUBLE PRECISION HM,HMUE,HEFF,ECRHMU,MODHMU
00297       DOUBLE PRECISION Q1,Q2,Q3
00298 !
00299       DOUBLE PRECISION CBID,FFW
00300       DOUBLE PRECISION PI,DEGRAD,RADDEG
00301 !
00302 !--> VARIABLES FOR CURRENT AND TETAP LOOP
00303       DOUBLE PRECISION ERREUR1,ERREURT
00304       DOUBLE PRECISION XMUL,XK,ZERO
00305       DOUBLE PRECISION TETA(NPTFR)   , ANGDIR(NPTFR)
00306 
00307       INTEGER ITERKN
00308 
00309 !
00310 !
00311 !-----------------------------------------------------------------------
00312 !
00313       PARAMETER( PI = 3.1415926535897932384626433D0 , DEGRAD=PI/180.D0 )
00314       PARAMETER( RADDEG = 180.D0 / PI )
00315 !
00316 !-----------------------------------------------------------------------
00317 !
00318       INTRINSIC ABS,MIN,MAX,LOG,COS,SIN
00319       DOUBLE PRECISION P_DMAX
00320       EXTERNAL P_DMAX
00321 
00322 
00323 !--------------------------
00324 !      LANGAUTO=.FALSE.
00325 !----------------------------------------------------------------------
00326 !
00327 ! INITIALISES MU AND FW: SET TO 0
00328 !         FOR THE FIRST ITERATION
00329 !
00330       ITERMU=0
00331       IF (LF.EQ.0) THEN
00332         CALL OS( 'X=C     ' , MU , SBID , SBID , 0.D0 )
00333         CALL OS( 'X=C     ' , FW , SBID , SBID , 0.D0 )
00334       ENDIF
00335 !
00336 !-----------------------------------------------------------------------
00337       ITERKN=0
00338       IF (COURANT) THEN
00339         ZERO=1E-10
00340 !       INITIALISATION OF WAVE VECTOR COMPONENTS X&Y : T5 , T6
00341         CALL OS('X=0     ',X=KANCX)
00342         CALL OS('X=0     ',X=KANCY)
00343       ENDIF
00344 !      WRITE(6,*) 'AVANT BOUCLE 98'
00345 
00346 !
00347 !-----------------------------------------------------------------------
00348 !
00349 !     =========================================
00350 !98 : DISSIPATION : ITERATIVE LOOP : ON THE VARIABLE MU R ON THE WAVE VECTOR
00351 98    CONTINUE
00352 !
00353 !   ITERATIVE LOOP ON TETAP (AUTOMATIC CALCULATION) AND WAVE VECTOR (WAVE-CURRENT)
00354       IF (    ((COURANT).AND.(ITERKN.GT.0))
00355      &    .OR.((LANGAUTO ).AND.(ITERKN.GT.0))  ) THEN
00356 
00357 !   ==> CURRENT :
00358         IF (COURANT) THEN
00359 !       COMPUTE WAVE VECTOR (FIRST ITERATION U=0, SO NO NEED TO DO THAT)
00360 !       ---------------------------------------------------------------
00361           DO I=1,NPOIN
00362             XK =K%R(I)
00363             CALL SOLVELAMBDA(XK,
00364      &                   UC%R(I),VC%R(I),KANCX%R(I),KANCY%R(I),H%R(I))
00365             K%R(I) =XK
00366 
00367             WR%R(I)=SQRT(GRAV*K%R(I)*TANH(K%R(I)*H%R(I)))
00368             C%R(I) =WR%R(I)/K%R(I)
00369             CG%R(I)=0.5D0*C%R(I)*
00370      &           (1.D0 + 2.D0*K%R(I)*H%R(I)/SINH(2.D0*K%R(I)*H%R(I)))
00371           ENDDO
00372 !
00373         ENDIF
00374 !       ------------------------------------------------------
00375 !       ==> AUTOMATIC ANGLES
00376         IF (LANGAUTO) THEN
00377 !       COMPUTE TETAP ON THE BOUNDARY (FIRST ITERATION TETAP GIVEN BY USER, SO NO NEED TO DO THAT)
00378 !       ---------------------------------------------------------------
00379 !        ==> RELAXATION ON TETAP IF NECESARY : TETAP=TETAP + C*(TETAP-TETAPM) (default : C=0)
00380           DO I=1,NPTFR
00381             TETAP%R(I)=TETAPM%R(I)+RELTP*(TETAP%R(I)-TETAPM%R(I))
00382             TETAP%R(I)=MIN(TETAP%R(I),90D0)
00383           ENDDO
00384 !         TETAP STORAGE IN TETAPM
00385 !         --------------------------
00386           CALL OS( 'X=Y     ' , X=TETAPM , Y=TETAP )
00387         ENDIF
00388 !
00389 !       ------------------------------------------------------
00390 !        ACTUALIZATION OF BOUNDARY CONDITIONS
00391 !           1/ CURRENT     : K HAS CHANGED
00392 !           2/ AUTO ANGLES : TETAP HAS CHANGED
00393 !       ----------
00394         CALL PHBOR
00395 !       ----------
00396 !      ------------------------------------------------------
00397 !     END OF FIRST STEP : NEXT STEP IS COMPUTING AM AND BM
00398       ENDIF
00399 !
00400 !     =========================================
00401 !                   MATRIX AM
00402 !     =========================================
00403 !      WRITE(6,*) 'MATRIX AM'
00404 !     ---------------------------
00405 !     DIFFUSION MATRIX FOR AM1
00406 !     ---------------------------
00407 !
00408       CALL OS( 'X=YZ    ' , T1 , C , CG , CBID )
00409       CALL MATRIX(AM1,'M=N     ','MATDIF          ',IELM,IELM,
00410      &            1.D0,S,S,S,T1,T1,S,MESH,MSK,MASKEL)
00411 !
00412 !-----------------------------------------------------------------------
00413 !
00414 ! PANCHANG, TO BE REVISITED: 7 IS GMRES
00415 !
00416 ! THE DIFFUSION MATRIX USED FOR PRECONDITIONING IS STORED
00417 ! IF THE METHOD IS THAT OF PANCHANG ET AL. (ISOLVE(1) =7)
00418 !
00419 !     IF (ISOLVE(1).EQ.7) THEN
00420 !
00421 !        CALL OM('M=CN    ',AM3,AM1,Z,1.D0/(RELAX*(2.D0-RELAX)),MESH)
00422 !
00423 !     ENDIF
00424 !
00425 !-----------------------------------------------------------------------
00426 !
00427 !     -----------------------
00428 !     MASS MATRIX FOR AM1
00429 !     -----------------------
00430 !     (WARNING : CAN'T USE CURRENT AND SECOND ORDER BOTTOM EFFECTS AT THE SAME TIME)
00431       IF (COURANT) THEN
00432         XMUL=1.D0
00433         CALL OS( 'X=YZ    ' , X=T1,Y=CG,Z=C)
00434         CALL OS( 'X=YZ    ' , X=T2 , Y=K , Z=K)
00435         CALL OS( 'X=XY    ' , X=T1 , Y=T2)
00436         CALL OS( 'X=C     ' , X=T2 , C=OMEGA**2)
00437         IF(ITERKN.EQ.0)THEN
00438           CALL OS( 'X=C     ' , X=T3 , C=OMEGA**2)
00439         ELSE
00440           CALL OS( 'X=YZ    ' , X=T3 , Y=WR , Z=WR )
00441         ENDIF
00442         CALL OS( 'X=Y+Z   ' , X=T1 , Y=T1 , Z=T2)
00443         CALL OS( 'X=Y-Z   ' , X=T1 , Y=T1 , Z=T3)
00444         IF (IPENTCO.GT.(0.5)) THEN
00445           WRITE(LU,*) 'IT IS NOT POSSIBLE TO USE '
00446           WRITE(LU,*) 'CURRENT + BOTTOM EFFECTS AT THE SAME TIME'
00447           WRITE(LU,*) '- FOR CURRENT ALONE, FIX IPENTO=0'
00448           WRITE(LU,*) '- FOR BOTTOM EFFECTS ALONE, DON T USE CURRENT'
00449           WRITE(LU,*) '-------------------------'
00450           WRITE(LU,*) 'THE CODE IS GOING TO STOP'
00451           WRITE(LU,*) '-------------------------'
00452           STOP
00453         ENDIF
00454       ELSE
00455         CALL OS( 'X=Y/Z   ' , X=T1,Y=CG,Z=C)
00456         XMUL=OMEGA**2
00457 !       SECOND ORDER BOTTOM EFFECTS ?
00458 !          (IPENTCO > 0 --> T1 = T1*(1+F) )
00459 !           0 : NO EFFECT /  1 : GRADIENT / 2 : CURVATURE /  3 : GRADIENT+CURVATURE
00460         IF ( (IPENTCO.GT.(0.5)).AND.(IPENTCO.LT.(3.5)) ) THEN
00461           CALL PENTCO(IPENTCO)
00462 !         WT USED AND TO BE CONSERVED : T3 = 1+F
00463 !         WT USED : T2 T4 T5 T6 T7 T9 T8 T11 T12
00464           CALL OS( 'X=YZ    ' , T1 , T1 , T3 , CBID )
00465         ENDIF
00466       ENDIF
00467 
00468       CALL MATRIX(AM2,'M=N     ','FMATMA          ', IELM , IELM ,
00469      &            XMUL , T1,S,S,S,S,S,MESH,MSK,MASKEL)
00470 !
00471 
00472 !     --------------------------------------------------
00473 !     COMPUTES DIFFUSION MATRIX - MASS MATRIX
00474 !     --------------------------------------------------
00475       CALL OM( 'M=M+CN  ' , AM1 , AM2 , C , -1.D0 , MESH )
00476 !
00477 
00478 !
00479       IF (COURANT) THEN
00480 !     --------------------------------------------------
00481 !     ADDS CURRENT - CONVECTION MATRIX
00482 !     --------------------------------------------------
00483 !
00484         CALL MATRIX(AM2,'M=N     ','MAUGUG          ', IELM , IELM ,
00485      &            1.D0 , S,S,S,UC,VC,S,MESH,MSK,MASKEL)
00486 !
00487         CALL OM( 'M=M+CN  ' , AM1 , AM2 , C , -1.D0 , MESH )
00488 !
00489       ENDIF
00490 !     --------------------------------
00491 !     ADDS THE BOUNDARY TERM TO AM1
00492 !     --------------------------------
00493 !
00494 !     AM1 et AM2 --> NON-SYMETRIQUES
00495       CALL OM( 'M=X(M)  ' , AM1 , AM1 , SBID , CBID , MESH )
00496       CALL OM( 'M=X(M)  ' , AM2 , AM2 , SBID , CBID , MESH )
00497 !
00498 !     ----------------------------
00499 !     BOUNDARY TERM: INCIDENT WAVE
00500 !     ----------------------------
00501 !
00502       IF (NPTFR .GT. 0) THEN
00503         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00504      &        -1.D0,BPHI1B,S,S,S,S,S,MESH,.TRUE.,MASK1)
00505         CALL OM( 'M=M+N   ' , AM1 , MBOR , T1 , CBID , MESH )
00506 !
00507 !
00508         IF(COURANT) THEN
00509           CALL MATRIX(MBOR,'M=N     ','MATFGUG         ',IELMB,IELMB,
00510      &        1.D0,MESH%XSGBOR,MESH%YSGBOR,S,UC,VC,S,
00511      &        MESH,.TRUE.,MASK1)
00512           CALL OM( 'M=M+N   ' , AM1 , MBOR , SBID , CBID , MESH )
00513         ENDIF
00514 !
00515       ENDIF
00516 !     ------------------------------
00517 !     BOUNDARY TERM: FREE EXIT
00518 !     ------------------------------
00519       IF (NPTFR .GT. 0) THEN
00520         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00521      &        -1.D0,BPHI2B,S,S,S,S,S,MESH,.TRUE.,MASK2)
00522         CALL OM( 'M=M+N   ' , AM1 , MBOR , T1 , CBID , MESH )
00523 !
00524         IF(COURANT) THEN
00525           CALL MATRIX(MBOR,'M=N     ','MATFGUG         ',IELMB,IELMB,
00526      &        1.D0,MESH%XSGBOR,MESH%YSGBOR,S,UC,VC,S,
00527      &        MESH,.TRUE.,MASK2)
00528           CALL OM( 'M=M+N   ' , AM1 , MBOR , SBID , CBID , MESH )
00529         ENDIF
00530       ENDIF
00531 !     ------------------------------
00532 !     BOUNDARY TERM : SOLID BOUNDARY
00533 !     ------------------------------
00534       IF (NPTFR .GT. 0) THEN
00535         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00536      &        -1.D0,BPHI3B,S,S,S,S,S,MESH,.TRUE.,MASK3)
00537         CALL OM( 'M=M+N   ' , AM1 , MBOR , T1 , CBID , MESH )
00538       END IF
00539 !     --------------------- ---------
00540 !     BOUNDARY TERM: IMPOSED WAVE
00541 !     ------------------------------
00542       IF (NPTFR .GT. 0) THEN
00543         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00544      &        -1.D0,BPHI4B,S,S,S,S,S,MESH,.TRUE.,MASK4)
00545         CALL OM( 'M=M+N   ' , AM1 , MBOR , T1 , CBID , MESH )
00546       END IF
00547 !
00548 !     ------------------------------
00549 !     BOUNDARY TERM: INCIDENT POTENTIAL
00550 !     ------------------------------
00551 !
00552       IF (NPTFR .GT. 0) THEN
00553         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00554      &        -1.D0,BPHI1B,S,S,S,S,S,MESH,.TRUE.,MASK5)
00555         CALL OM( 'M=M+N   ' , AM1 , MBOR , T1 , CBID , MESH )
00556 !
00557         IF(COURANT) THEN
00558           CALL MATRIX(MBOR,'M=N     ','MATFGUG         ',IELMB,IELMB,
00559      &        1.D0,MESH%XSGBOR,MESH%YSGBOR,S,UC,VC,S,
00560      &        MESH,.TRUE.,MASK5)
00561           CALL OM( 'M=M+N   ' , AM1 , MBOR , SBID , CBID , MESH )
00562         ENDIF
00563 !
00564       ENDIF
00565 !
00566 !     =========================================
00567 !                   SECOND MEMBERS
00568 !     =========================================
00569 !     WRITE(6,*) 'CV1 et CV2'
00570 !     ---------------------
00571 !     SECOND MEMBERS : CV1
00572 !     ---------------------
00573 !
00574       CALL OS( 'X=C     ' , CV1, SBID , SBID , 0.D0 )
00575 !     ------------------------------
00576 !     BOUNDARY TERM: INCIDENT WAVE
00577 !     ------------------------------
00578 !       --- CALCUL DE i COS(TETAP) GAMMA
00579       CALL OS( 'X=CY    ' , T1,TETAP,SBID,DEGRAD)
00580       CALL OS( 'X=COS(Y)' , T2,T1,SBID,0.D0)
00581       CALL OS( 'X=YZ    ' , T3,CPHI1B,T2,0.D0)
00582       CALL OS( 'X=C     ' , T1, SBID , SBID , 0.D0 )
00583 
00584       IF (NPTFR .GT. 0) THEN
00585         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00586      &           -1.D0,T3,S,S,S,S,S,MESH,.TRUE.,MASK1)
00587       END IF
00588       CALL OSDB( 'X=X+Y   ' , CV1 , T1 , SBID , CBID , MESH )
00589 
00590 !     --- CALCUL DE GRAD(Gamma).n : REEL
00591       CALL OS( 'X=Y     ' , T2,CGRX1B,SBID,0.D0)
00592       CALL OS( 'X=Y     ' , T3,CGRY1B,SBID,0.D0)
00593       CALL OS( 'X=C     ' , T1, SBID , SBID , 0.D0 )
00594       CALL OS( 'X=C     ' , T4, SBID , SBID , 1.D0 )
00595       IF (NPTFR .GT. 0) THEN
00596         CALL VECTOR(T1,'=','FLUBDF          ',IELMB,
00597      &           1.D0,T4,S,S,T2,T3,S,MESH,.TRUE.,MASK1)
00598       END IF
00599       CALL OSDB( 'X=X+Y   ' , CV1 , T1 , SBID , CBID , MESH )
00600 !     ---------------------------------
00601 !     BOUNDARY TERM: INCIDENT POTENTIAL
00602 !     ---------------------------------
00603 !     --- CALCUL DE i COS(TETAP) GAMMA
00604       CALL OS( 'X=CY    ' , T1,TETAP,SBID,DEGRAD)
00605       CALL OS( 'X=COS(Y)' , T2,T1,SBID,0.D0)
00606       CALL OS( 'X=YZ    ' , T3,CPHI1B,T2,0.D0)
00607       CALL OS( 'X=C     ' , T1, SBID , SBID , 0.D0 )
00608       IF (NPTFR .GT. 0) THEN
00609         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00610      &           -1.D0,T3,S,S,S,S,S,MESH,.TRUE.,MASK5)
00611       END IF
00612       CALL OSDB( 'X=X+Y   ' , CV1 , T1 , SBID , CBID , MESH )
00613 
00614 !     --- CALCUL DE GRAD(Gamma).n : REEL
00615       CALL OS( 'X=Y     ' , T2,CGRX1B,SBID,0.D0)
00616       CALL OS( 'X=Y     ' , T3,CGRY1B,SBID,0.D0)
00617       CALL OS( 'X=C     ' , T1, SBID , SBID , 0.D0 )
00618       CALL OS( 'X=C     ' , T4, SBID , SBID , 1.D0 )
00619       IF (NPTFR .GT. 0) THEN
00620         CALL VECTOR(T1,'=','FLUBDF          ',IELMB,
00621      &           1.D0,T4,S,S,T2,T3,S,MESH,.TRUE.,MASK5)
00622       END IF
00623       CALL OSDB( 'X=X+Y   ' , CV1 , T1 , SBID , CBID , MESH )
00624 !     ------------------------------
00625 !     BOUNDARY TERM: FREE EXIT
00626 !     ------------------------------
00627       IF (NPTFR .GT. 0) THEN
00628         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00629      &           1.D0,CPHI2B,S,S,S,S,S,MESH,.TRUE.,MASK2)
00630         CALL OSDB( 'X=X+Y   ' , CV1 , T1 , SBID , CBID , MESH )
00631       END IF
00632 !     ------------------------------
00633 !     BOUNDARY TERM: SOLID BOUNDARY
00634 !     ------------------------------
00635       IF (NPTFR .GT. 0) THEN
00636         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00637      &          1.D0,CPHI3B,S,S,S,S,S,MESH,.TRUE.,MASK3)
00638         CALL OSDB( 'X=X+Y   ' , CV1 , T1 , SBID , CBID , MESH )
00639       END IF
00640 !     ------------------------------
00641 !     BOUNDARY TERM: IMPOSED WAVE
00642 !     ------------------------------
00643       IF (NPTFR .GT. 0) THEN
00644         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00645      &           1.D0,CPHI4B,S,S,S,S,S,MESH,.TRUE.,MASK4)
00646         CALL OSDB( 'X=X+Y   ' , CV1 , T1 , SBID , CBID , MESH )
00647       END IF
00648 !
00649 !CP
00650 !         IF (NCSIZE.GT.1) THEN
00651 !           CALL PARCOM(CV1,2,MESH)
00652 !         ENDIF
00653 !CP
00654 !     ---------------------
00655 !     SECOND MEMBERS : CV2
00656 !     ---------------------
00657 !
00658       CALL OS( 'X=C     ' , CV2, SBID , SBID , 0.D0 )
00659 !     ------------------------------
00660 !     BOUNDARY TERM: INCIDENT WAVE
00661 !     ------------------------------
00662 !     --- CALCUL DE i COS(TETAP) GAMMA : IMAGINAIRE
00663       CALL OS( 'X=CY    ' , T1,TETAP,SBID,DEGRAD)
00664       CALL OS( 'X=COS(Y)' , T2,T1,SBID,0.D0)
00665       CALL OS( 'X=YZ    ' , T3,DPHI1B,T2,0.D0)
00666       CALL OS( 'X=C     ' , T1, SBID , SBID , 0.D0 )
00667       IF (NPTFR .GT. 0) THEN
00668         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00669      &           -1.D0,T3,S,S,S,S,S,MESH,.TRUE.,MASK1)
00670       END IF
00671       CALL OSDB( 'X=X+Y   ' , CV2 , T1 , SBID , CBID , MESH )
00672 
00673 !     --- CALCUL DE GRAD(Gamma).n : IMAGINAIRE
00674       CALL OS( 'X=Y     ' , T2,DGRX1B,SBID,0.D0)
00675       CALL OS( 'X=Y     ' , T3,DGRY1B,SBID,0.D0)
00676       CALL OS( 'X=C     ' , T1, SBID , SBID , 0.D0 )
00677       CALL OS( 'X=C     ' , T4, SBID , SBID , 1.D0 )
00678 
00679       IF (NPTFR .GT. 0) THEN
00680         CALL VECTOR(T1,'=','FLUBDF          ',IELMB,
00681      &           1.D0,T4,S,S,T2,T3,S,MESH,.TRUE.,MASK1)
00682       END IF
00683       CALL OSDB( 'X=X+Y   ' , CV2 , T1 , SBID , CBID , MESH )
00684 !     ---------------------------------
00685 !     BOUNDARY TERM: INCIDENT POTENTIAL
00686 !     ---------------------------------
00687 !     --- CALCUL DE i COS(TETAP) GAMMA : IMAGINAIRE
00688       CALL OS( 'X=CY    ' , T1,TETAP,SBID,DEGRAD)
00689       CALL OS( 'X=COS(Y)' , T2,T1,SBID,0.D0)
00690       CALL OS( 'X=YZ    ' , T3,DPHI1B,T2,0.D0)
00691       CALL OS( 'X=C     ' , T1, SBID , SBID , 0.D0 )
00692       IF (NPTFR .GT. 0) THEN
00693         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00694      &           -1.D0,T3,S,S,S,S,S,MESH,.TRUE.,MASK5)
00695       END IF
00696       CALL OSDB( 'X=X+Y   ' , CV2 , T1 , SBID , CBID , MESH )
00697 
00698 !     --- CALCUL DE GRAD(Gamma).n : IMAGINAIRE
00699       CALL OS( 'X=Y     ' , T2,DGRX1B,SBID,0.D0)
00700       CALL OS( 'X=Y     ' , T3,DGRY1B,SBID,0.D0)
00701       CALL OS( 'X=C     ' , T1, SBID , SBID , 0.D0 )
00702       CALL OS( 'X=C     ' , T4, SBID , SBID , 1.D0 )
00703 
00704       IF (NPTFR .GT. 0) THEN
00705         CALL VECTOR(T1,'=','FLUBDF          ',IELMB,
00706      &           1.D0,T4,S,S,T2,T3,S,MESH,.TRUE.,MASK5)
00707       END IF
00708       CALL OSDB( 'X=X+Y   ' , CV2 , T1 , SBID , CBID , MESH )
00709 !
00710 !    ------------------------------
00711 !    BOUNDARY TERM: FREE EXIT
00712 !    ------------------------------
00713       IF (NPTFR .GT. 0) THEN
00714         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00715      &           1.D0,DPHI2B,S,S,S,S,S,MESH,.TRUE.,MASK2)
00716       END IF
00717       CALL OSDB( 'X=X+Y   ' , CV2 , T1 , SBID , CBID , MESH )
00718 !     ------------------------------
00719 !     BOUNDARY TERM: SOLID BOUNDARY
00720 !     -----------------------------
00721       IF (NPTFR .GT. 0) THEN
00722         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00723      &           1.D0,DPHI3B,S,S,S,S,S,MESH,.TRUE.,MASK3)
00724         CALL OSDB( 'X=X+Y   ' , CV2 , T1 , SBID , CBID , MESH )
00725       END IF
00726 !     ------------------------------
00727 !     BOUNDARY TERM: IMPOSED WAVE
00728 !     ------------------------------
00729       IF (NPTFR .GT. 0) THEN
00730         CALL VECTOR(T1,'=','MASVEC          ',IELMB,
00731      &           1.D0,DPHI4B,S,S,S,S,S,MESH,.TRUE.,MASK4)
00732         CALL OSDB( 'X=X+Y   ' , CV2 , T1 , SBID , CBID , MESH )
00733       END IF
00734 !CP
00735 !          IF (NCSIZE.GT.1) THEN
00736 !           CALL PARCOM(CV2,2,MESH)
00737 !         ENDIF
00738 !CP
00739 !     =========================================
00740 !                   MATRIX BM
00741 !     =========================================
00742 !      WRITE(6,*) 'MATRIX BM'
00743 !
00744 !     ----------------------------------------------------------
00745 !     COMPUTES THE MATRIX BM1 FOR THE MU VALUES SPECIFIED
00746 !     FOR THE ITERATION 'ITERMU'
00747 !     ----------------------------------------------------------
00748 !
00749       CALL OS( 'X=YZ    ' , T1 , C  , CG , CBID )
00750       CALL OS( 'X=YZ    ' , T2 , K  , MU , CBID )
00751       CALL OS( 'X=YZ    ' , T1 , T1 , T2 , CBID )
00752       CALL MATRIX(BM1,'M=N     ','FMATMA          ', IELM , IELM ,
00753      &            1.D0 , T1,S,S,S,S,S,MESH,MSK,MASKEL)
00754 !
00755       IF ((COURANT).AND.(ITERKN.GT.0)) THEN
00756 !     ----------------------------------------------------------
00757 !     ADD TERMS TO BM1 FOR THE CURRENT 2*OMEGA.U TERM
00758 !     ----------------------------------------------------------
00759 !          ON DESYMETRISE BM1
00760         CALL OM( 'M=X(M)  ' , BM1 , BM1 , SBID , CBID , MESH )
00761 !
00762         CALL MATRIX(BM2,'M=N     ','MATVGR          ',IELM ,IELM ,
00763      &            2D0*OMEGA , S,S,S,UC,VC,S,
00764      &            MESH,MSK,MASKEL)
00765 !
00766         CALL OM( 'M=M+N   ' , BM1 , BM2 , T1 , CBID , MESH )
00767 !
00768 !     ----------------------------------------------------------
00769 !     ADD TERMS TO THE MATRIX BM1 FOR THE CURRENT : DIV(U) TERM
00770 !     ----------------------------------------------------------
00771 !
00772         CALL VECTOR(T2 , '=' , 'MASBAS          ' , IELM ,
00773      &            1.D0 , S , S , S , S , S , S ,
00774      &            MESH , MSK  , MASKEL )
00775 !
00776         CALL VECTOR(T13 , '=' , 'GRADF          X' , IELM ,
00777      &            1.D0 , UC , S , S , S , S , S ,
00778      &            MESH , MSK , MASKEL)
00779 !
00780         CALL VECTOR(T14 , '=' , 'GRADF          Y' , IELM ,
00781      &            1.D0 , VC , S , S , S , S , S ,
00782      &            MESH , MSK , MASKEL)
00783 !
00784         CALL OS( 'X=Y+Z   ' , X=T15 , Y=T14 , Z=T13 )
00785         CALL OS( 'X=Y/Z   ' , X=T16 , Y=T15 , Z=T2 )
00786 !
00787         CALL MATRIX(BM2,'M=N     ','FMATMA          ',IELM ,IELM ,
00788      &            OMEGA , T16,S,S,S,S,S,
00789      &            MESH,MSK,MASKEL)
00790 !
00791         CALL OM( 'M=M+N   ' , BM1 , BM2 , SBID , CBID , MESH )
00792       ENDIF
00793 
00794 
00795 !     -------------------------------------------
00796 !     ADDS THE BOUNDARY TERM TO BM1
00797 !     -------------------------------------------
00798 !
00799       IF (NPTFR .GT. 0) THEN
00800 !        ------------------------------
00801 !        BOUNDARY TERM: INCIDENT WAVE
00802 !        ------------------------------
00803         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00804      &        -1.D0,APHI1B,S,S,S,S,S,MESH,.TRUE.,MASK1)
00805 !         WRITE (*,*)'MBOR = ', MBOR%TYPEXT
00806         CALL OM( 'M=M+N   ' , BM1 , MBOR , T1 , CBID , MESH )
00807       END IF
00808 !        ------------------------------
00809 !        BOUNDARY TERM: FREE EXIT
00810 !        ------------------------------
00811       IF (NPTFR .GT. 0) THEN
00812         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00813      &           -1.D0,APHI2B,S,S,S,S,S,MESH,.TRUE.,MASK2)
00814         CALL OM( 'M=M+N   ' , BM1 , MBOR , T1 , CBID , MESH )
00815       END IF
00816 !        ------------------------------
00817 !        BOUNDARY TERM: SOLID BOUNDARY
00818 !        ------------------------------
00819       IF (NPTFR .GT. 0) THEN
00820         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00821      &        -1.D0,APHI3B,S,S,S,S,S,MESH,.TRUE.,MASK3)
00822         CALL OM( 'M=M+N   ' , BM1 , MBOR , T1 , CBID , MESH )
00823       END IF
00824 !        ------------------------------
00825 !        BOUNDARY TERM: IMPOSED WAVE
00826 !        ------------------------------
00827       IF (NPTFR .GT. 0) THEN
00828         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00829      &           -1.D0,APHI4B,S,S,S,S,S,MESH,.TRUE.,MASK4)
00830         CALL OM( 'M=M+N   ' , BM1 , MBOR , T1 , CBID , MESH )
00831       END IF
00832 !        ------------------------------
00833 !        BOUNDARY TERM: INCIDENT POTENTIAL
00834 !        ------------------------------
00835       IF (NPTFR .GT. 0) THEN
00836         CALL MATRIX(MBOR,'M=N     ','FMATMA          ',IELMB,IELMB,
00837      &        -1.D0,APHI1B,S,S,S,S,S,MESH,.TRUE.,MASK5)
00838         CALL OM( 'M=M+N   ' , BM1 , MBOR , T1 , CBID , MESH )
00839       END IF
00840 
00841 !     ---------
00842 !     AM2 = AM1
00843 !     ---------
00844 !
00845       CALL OM( 'M=N     ' , AM2 , AM1 , SBID , CBID , MESH )
00846 !
00847 !     --------------------------
00848 !     BM1 BECOMES NONSYMMETRICAL
00849 !     --------------------------
00850 !
00851 !     SI ON N'A PAS DE COURANT OU SI ON NE L'A PAS ENCORE PRIS EN COMPTE
00852 !
00853 !      WRITE(6,*) 'MATRIX BM non sym'
00854       IF ((.NOT.COURANT).OR.ITERKN.EQ.0) THEN
00855         CALL OM( 'M=X(M)  ' , BM1 , BM1 , SBID , CBID , MESH )
00856       ENDIF
00857 
00858 !
00859 !     ----------------------------
00860 !     TRIES MASS-LUMPING OF BM1
00861 !     ----------------------------
00862 !
00863 !     MASLU = 1.D0
00864 !     CALL LUMP(T1,BM1,MESH,XMESH,MASLU,MSK,MASKEL)
00865 !     CALL OM( 'M=CN    ' , BM1 , BM1 , T1 , 1.D0-MASLU , MESH )
00866 !     CALL OM( 'M=M+D   ' , BM1 , BM1 , T1 , C          , MESH )
00867 !
00868 !     ----------
00869 !     BM2 = -BM1
00870 !     ----------
00871 !
00872       CALL OM( 'M=CN    ' , BM2 , BM1 , C , -1.D0 , MESH )
00873 !
00874 !     =======================================
00875 !
00876 !     TAKES INTO ACCOUNT DIRICHLET POINTS
00877 !
00878 !     =======================================
00879 !
00880       IF ((.NOT.ALEMON).AND.(.NOT.ALEMUL)) THEN
00881         IF (DEFERL .OR. FROTTE) THEN
00882           IF (LNG.EQ.1) WRITE(LU,220) ITERMU+1
00883           IF (LNG.EQ.2) WRITE(LU,221) ITERMU+1
00884         ENDIF
00885       ENDIF
00886  220  FORMAT(/,1X,'SOUS-ITERATION NUMERO :',1X,I3,/)
00887  221  FORMAT(/,1X,'SUB-ITERATION NUMBER :',1X,I3,/)
00888 !
00889       CALL DIRICH(UNK,MAT,RHS,PHIB,LIDIR%I,TB,MESH,KENT,MSK,MASKEL)
00890 !
00891 !     ===============================================================
00892 !
00893 !     INHIBITS POSSIBLE DIAGONAL PRECONDITIONING
00894 !     IF AN ELEMENT OF DAM1 IS NEGATIVE OR NULL
00895 !
00896 !     ===============================================================
00897 !
00898       CALL CNTPRE(AM1%D%R,NPOIN,SLVART%PRECON,SLVART%PRECON)
00899 !      IF (LNG.EQ.1) WRITE(LU,230) SLVART%PRECON
00900 !      IF (LNG.EQ.2) WRITE(LU,231) SLVART%PRECON
00901 ! 230  FORMAT(/,1X,'PRECONDITIONNEMENT APRES CONTROLE :',1X,I3)
00902 ! 231  FORMAT(/,1X,'PRECONDITIONNING AFTER CONTROL :',1X,I3)
00903 !
00904 !     ==========================================================
00905 !
00906 !     PRECONDITIONING BLOCK-DIAGONAL:
00907 !                 THE MATRICES BECOME NONSYMMETRICAL.
00908 !
00909 !     ==========================================================
00910 !
00911       IF (3*(SLVART%PRECON/3).EQ.SLVART%PRECON) THEN
00912         CALL OM( 'M=X(M)  ' , AM1 , AM1 , SBID , CBID , MESH )
00913         CALL OM( 'M=X(M)  ' , AM2 , AM2 , SBID , CBID , MESH )
00914       ENDIF
00915 !
00916 !     ==============================
00917 !
00918 !     SOLVES THE LINEAR SYSTEM
00919 !
00920 !     ==============================
00921 !
00922 !     ----------------------------
00923 !     INITIALISES THE UNKNOWN
00924 !     ----------------------------
00925 !
00926       IF(ITERMU.EQ.0.AND.LF.EQ.0) THEN
00927         CALL LUMP(T1,AM1,MESH,1.D0)
00928         CALL OS( 'X=Y/Z   ' , PHIR , CV1 , T1 , CBID )
00929         CALL LUMP(T1,AM2,MESH,1.D0)
00930         CALL OS( 'X=Y/Z   ' , PHII , CV2 , T1 , CBID )
00931       ENDIF
00932 !
00933       IF (LNG.EQ.1) WRITE(LU,240)
00934       IF (LNG.EQ.2) WRITE(LU,241)
00935  240  FORMAT(/,1X,'RESOLUTION DU SYSTEME LINEAIRE (SOLVE)',/)
00936  241  FORMAT(/,1X,'LINEAR SYSTEM SOLVING (SOLVE)',/)
00937 !
00938       CALL SOLVE(UNK,MAT,RHS,TB,SLVART,INFOGR,MESH,AM3)
00939 !
00940 !     ============================================================
00941 !     DIRECTION LOOP
00942 !      - WAVE-CURRENT :checks convergence on the wave vector
00943 !      - AUTO ANGLES  :checks convergence on TETAP
00944 !     ============================================================
00945       IF (COURANT.OR.LANGAUTO) THEN
00946 !      ------------------------------------------------------
00947 !      COMPUTE WAVE INCIDENCE USING SPEED AT THE FREE SURFACE
00948 !
00949 !       -----------
00950         CALL CALDIR()
00951 !       -----------
00952 !          --> PHIR,PHII
00953 !          --  T1,T2,T3,T4
00954 !         <--  INCI
00955 !        -- DIRECTION OF VECTOR K : INCI
00956         CALL OS( 'X=COS(Y)' , T5,INCI,SBID,0.D0)
00957         CALL OS( 'X=SIN(Y)' , T6,INCI,SBID,0.D0)
00958 !       T5 = K  COS(INCIDENCE)
00959         CALL OS( 'X=XY    ' , X=T5 , Y=K)
00960 !       T6 = K  SIN(INCIDENCE)
00961         CALL OS( 'X=XY    ' , X=T6 , Y=K)
00962 !       ------------------------------------------------------
00963 !
00964 !       Error Initialisation
00965         ERREUR1=0.D0
00966         ERREURT=0.D0
00967 !     --------------------------------------------------
00968 !     CONVERGENCE CRITERION FOR WAVE-CURRENT INTERACTION
00969         IF (COURANT) THEN
00970 !       MAX ERROR CALCULATION : NORM( Kn - Kn-1 )/NORM(Kn) < EPSDIR
00971 !       ----------------------
00972           IF (ITERKN.GT.0) THEN
00973             DO I=1,NPOIN
00974               ERREUR1=MAX(ERREUR1,
00975      &          SQRT((KANCX%R(I)-T5%R(I))**2+(KANCY%R(I)-T6%R(I))**2)/
00976      &          SQRT(T5%R(I)**2+T6%R(I)**2)
00977      &               )
00978             ENDDO
00979             WRITE(LU,*) '--------------------------------------------'
00980             WRITE(LU,*) 'WAVE-CURRENT : DIFF. BETWEEN 2 ITER. =',
00981      &                   ERREUR1
00982             WRITE(LU,*) 'LOOP FOR WAVE-CURRENT : TOLERANCE    =',
00983      &                   EPSDIR
00984           ELSE
00985             WRITE(LU,*) 'INITIAL LOOP FOR WAVE-CURRENT COMPLETED'
00986           ENDIF
00987           WRITE(LU,*) '----------------------------------------------'
00988 !       OLD WAVE VECTOR STORAGE
00989 !       -----------------------
00990           CALL OS( 'X=Y     ' , X=KANCX , Y=T5 )
00991           CALL OS( 'X=Y     ' , X=KANCY , Y=T6 )
00992         ENDIF
00993 !     --------------------------------------------------
00994 !
00995 !     -----------------------------------------------------
00996 !     CONVERGENCE CRITERION FOR TETAP AUTOMATIC CALCULATION
00997         IF (LANGAUTO) THEN
00998 !       STORAGE OF INCIDENCE ANGLE ON THE BOUNDARY IN A TABLE
00999           DO I=1,NPTFR
01000             IG       = MESH%NBOR%I(I)
01001             ANGDIR(I)=INCI%R(IG)
01002           ENDDO
01003 !         TETAP COMPUTATION
01004           CALL CALTETAP(TETA,MESH%XNEBOR%R,MESH%YNEBOR%R,
01005      &                     MESH%XSGBOR%R,MESH%YSGBOR%R,ANGDIR,NPTFR)
01006 !       MAX ERROR CALCULATION : MAX(cos(TETAPnew) - cos(TETAPold)) < EPSTP
01007 !       ----------------------
01008           IF (ITERKN.GT.0) THEN
01009             DO I=1,NPTFR
01010               TETAP%R(I)=TETA(I)
01011 !             ADD FACTOR 1-R%P EXP(iALFA) ?  ADD PHI?
01012               ERREURT=MAX(ERREURT,
01013      &             ABS(COS(TETAP%R(I)*DEGRAD)-COS(TETAPM%R(I)*DEGRAD))
01014      &                )
01015             ENDDO
01016             WRITE(LU,*) '-------------------------------------------'
01017             WRITE(LU,*) 'AUTO-ANGLES : DIFF. BETWEEN 2 ITER. =',
01018      &                   ERREURT
01019             WRITE(LU,*) 'LOOP FOR AUTO-ANGLE  : TOLERANCE    =',
01020      &                   EPSTP
01021           ELSE
01022             DO I=1,NPTFR
01023               TETAPM%R(I)=TETAP%R(I)
01024               TETAP%R(I) =TETA(I)
01025             ENDDO
01026             WRITE(LU,*) 'INITIAL LOOP FOR AUTOMATIC ANGLES COMPLETED'
01027           ENDIF
01028           WRITE(LU,*) '--------------------------------------------'
01029         ENDIF
01030 !     -----------------------------------------------------
01031 !
01032 !
01033 !      MAX ERROR FOR N PROC
01034         IF (NCSIZE.GT.1) THEN
01035           ERREURT = P_DMAX(ERREURT)
01036           ERREUR1 = P_DMAX(ERREUR1)
01037         END IF
01038 !
01039 !      ----------------------------------------------------
01040 !      CHECK CONVERGENCE FOR DIRECTION LOOP
01041         IF ( (ERREUR1.GT.EPSDIR).OR.(ERREURT.GT.EPSTP)
01042      &                         .OR.(ITERKN.EQ.0)     ) THEN
01043 !         NEW ITERATION
01044           ITERKN = ITERKN + 1
01045           IF (ITERKN.LE.NITTP) THEN
01046             GOTO 98
01047           ELSE
01048             IF (LNG.EQ.1) WRITE(LU,100) ITERKN
01049             IF (LNG.EQ.2) WRITE(LU,101) ITERKN
01050           ENDIF
01051         ENDIF
01052         IF (LNG.EQ.1) WRITE(LU,202) ITERKN
01053         IF (LNG.EQ.2) WRITE(LU,203) ITERKN
01054 !       REMISE A 1 DU NOMBRE d'ITERATION SUR LE COURRANT ET LA DIRECTION
01055         ITERKN=1
01056 !       =================================================
01057 !       END OF THE LOOP ON DIRECTIONS AND WAVE NUMBER
01058 !       =================================================
01059       ENDIF
01060 !    ----------------------------------------------------
01061 !
01062  100  FORMAT(/,1X,'BERKHO (ARTEMIS): NOMBRE DE SOUS-ITERATIONS',
01063      & 1X,'MAXIMUM ATTEINT POUR LE COURANT OU TETAP:',1X,I3)
01064  101  FORMAT(/,1X,'BERKHO (ARTEMIS): YOU REACHED THE MAXIMUM',
01065      & 1X,'NUMBER OF SUB-ITERATIONS FOR CURRENT OR TETAP :)',1X,I3)
01066 
01067  202  FORMAT(/,1X,'NOMBRE DE SOUS-ITERATIONS DIRECTION / COURANT :',
01068      &   1X,I3)
01069  203  FORMAT(/,1X,'NUMBER OF SUB-ITERATIONS DIRECTION / CURRENT :',
01070      &   1X,I3)
01071 
01072 !     ============================================================
01073 !
01074 !     COMPUTES THE TOTAL DISSIPATION COEFFICIENT MU_DEFERL + MU_FROTTE
01075 !     FOR REGULAR WAVES
01076 !     DISSIPATION FOR IRREGULAR WAVE IS LOOKED AT INTO ARTEMIS.F)
01077 !     ============================================================
01078 !
01079       IF (.NOT. ALEMON .AND. .NOT. ALEMUL) THEN
01080         IF (DEFERL .OR. FROTTE) THEN
01081           ECRHMU=0D0
01082 !         COMPUTES DISSIPATION COEFFICIENT MU2
01083           CALL CALCMU(ITERMU)
01084 !              WORK TABLE USED                      : T1,T4
01085 !              WORK TABLE USED AND TO BE CONSERVED  : T3 => QB
01086 !
01087 !         USE RELAXATION METHOD FOR DISSPATION COEFFICIENT MU
01088           CALL RELAXMU(ECRHMU,MODHMU,ITERMU)
01089 !              WORK TABLE USED                      : NONE
01090 
01091 !         ----------------------------------------------------
01092 !         CHECKS CONVERGENCE ON THE DISSIPATION ITERATIVE LOOP
01093 !         ----------------------------------------------------
01094           WRITE(LU,*) ' '
01095           WRITE(LU,*) '----------------------------------------------- '
01096           IF (ECRHMU.GT.EPSDIS*MODHMU) GOTO 98
01097 !
01098 !         QB STORAGE
01099 !         ----------
01100           CALL OS( 'X=Y     ', QB,T3,SBID,CBID)
01101 !
01102           IF (LNG.EQ.1) WRITE(LU,200) ITERMU
01103           IF (LNG.EQ.2) WRITE(LU,201) ITERMU
01104  200      FORMAT(/,1X,'NOMBRE DE SOUS-ITERATIONS POUR LA DISSIPATION:',
01105      &    1X,I3)
01106  201      FORMAT(/,1X,'NUMBER OF SUB-ITERATIONS FOR DISSIPATION:',
01107      &    1X,I3)
01108 !
01109         ENDIF
01110       ENDIF
01111 !
01112 !
01113 !-----------------------------------------------------------------------
01114 !
01115       RETURN
01116       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0