intemp.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\intemp.f
00002 !
00046                        SUBROUTINE INTEMP
00047 !                      *****************
00048 !
00049      &(W,FLUX,FLUX_OLD,AIRS,DT,NPOIN,ZF,CF,EPS,KFROT,SMH,
00050      & HN,QU,QV,LT,GAMMA)
00051 !
00052 !***********************************************************************
00053 ! TELEMAC2D   V6P1                                         03/15/2011
00054 !***********************************************************************
00055 !
00056 !           U_(N+1)=U_N + DT*( (1-GAMMA)ACC_N +GAMMA*ACC_(N+1))
00057 !       ACC: IS THE ACCELERATION (FLUX BALANCE FOR FV)
00058 !       FOR GAMMA=0.5 THE SCHEME IS SECOND ORDER ACCURATE
00059 !       FOR GAMMA=1.0 THE SCHEME IS EULER EXPLICIT (FIRST ORDER)
00060 !+
00061 !
00062 !
00063 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00064 !|  W             |<--|  (H,HU,HV)
00065 !|  FLUX          |-->|  FLUX AT TN+1
00066 !|  FLUX_OLD      |-->|  FLUX AT TN
00067 !|  AIRS          |-->|  CELL AREAS
00068 !|  DT            |-->|  TIME STEP
00069 !|  NPOIN         |-->|  TOTAL NUMBER OF NODES
00070 !|  ZF            |-->|  BATHYMETRIES
00071 !|  CF            |-->|  FRICTION COEFFICIENTS
00072 !|  EPS           |-->|  TOLERANCE FOR WATER DEPTH
00073 !|  KFROT         |-->|  LOGICAL! FRICTION OF NO FRICTION
00074 !|  SMH           |-->|  MASS SOURCE
00075 !|  HN,QU;QV      |-->|  H, HU AND HV AT TN
00076 !|  LT            |-->|  CURRENT TIME ITERATION
00077 !|  GAMMA         |-->|  NEWMARK PARAMETER (SEE BELOW)
00078 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00079 !
00080       IMPLICIT NONE
00081 !
00082       INTEGER LNG,LU
00083       COMMON/INFO/LNG,LU
00084 !
00085 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00086 !
00087       INTEGER, INTENT(IN)             :: NPOIN,KFROT,LT
00088       DOUBLE PRECISION, INTENT(INOUT) :: W(3,NPOIN)
00089       DOUBLE PRECISION, INTENT(IN)    :: FLUX(NPOIN,3),DT,EPS
00090       DOUBLE PRECISION, INTENT(IN)    :: AIRS(NPOIN),ZF(NPOIN)
00091       DOUBLE PRECISION, INTENT(IN)    :: FLUX_OLD(NPOIN,3)
00092       DOUBLE PRECISION, INTENT(IN)    :: CF(NPOIN),SMH(NPOIN),GAMMA
00093       DOUBLE PRECISION, INTENT(IN)    :: HN(NPOIN),QU(NPOIN),QV(NPOIN)
00094 !
00095 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00096 !
00097       INTEGER I
00098 !
00099       DOUBLE PRECISION DELTA,KF,ST2D,ALPHAF,UNMGAMMA,FACT
00100 !
00101 !=======================
00102 !---- NEWMARK SCHEME
00103 !=======================
00104 !
00105       UNMGAMMA = 1.D0-GAMMA
00106 !
00107 !     - FOR GAMMA=0.5, THIS CHOICE GIVES ORDER 2 ACCURACY AND
00108 !     THE SCHEME IS UNCONDITIALLY STABLE
00109 !     - FOR USER WHO PREFERS (EULER) EXPLICIT SCHEME,
00110 !     YOU HAVE TO PUT GAMMA=1
00111 !
00112       DO I=1,NPOIN
00113         FACT=DT/AIRS(I)
00114 !---    FIRST TIME STEP
00115         IF(LT.EQ.1)THEN
00116           W(1,I) = HN(I) - FACT * ( FLUX(I,1)-SMH(I) )
00117           W(2,I) = QU(I) - FACT * FLUX (I,2)
00118           W(3,I) = QV(I) - FACT * FLUX (I,3)
00119         ELSE
00120           W(1,I) = HN(I) - FACT * (UNMGAMMA*FLUX_OLD(I,1) +
00121      &                             GAMMA*FLUX(I,1)-SMH(I))
00122           W(2,I) = QU(I) - FACT * (UNMGAMMA*FLUX_OLD(I,2) +
00123      &                             GAMMA*FLUX(I,2))
00124           W(3,I) = QV(I) - FACT * (UNMGAMMA*FLUX_OLD(I,3) +
00125      &                             GAMMA*FLUX(I,3))
00126         ENDIF
00127       ENDDO
00128 !
00129 !     TO TAKE INTO ACCOUNT FRICTION
00130 !     *****************************
00131       IF (KFROT.NE.0) THEN
00132 !
00133         DO I = 1,NPOIN
00134 !
00135           IF((W(1,I).GT.EPS/10.D0).AND.(CF(I).GT.1.D-12)) THEN
00136             ST2D = CF(I)
00137             KF = 9.81D0*DT*DSQRT(W(2,I)**2+W(3,I)**2)/
00138      &           (ST2D*ST2D*W(1,I)**(7.D0/3.D0))
00139             IF(KF.GT.1.D-6) THEN
00140               DELTA = (1.D0+4.D0*KF)
00141               ALPHAF = (-1.D0+SQRT(DELTA))/(2*KF)
00142             ELSE
00143               ALPHAF = 1.D0 - KF
00144             ENDIF
00145             W(2,I) = ALPHAF * W(2,I)
00146             W(3,I) = ALPHAF * W(3,I)
00147           ENDIF
00148 !
00149         ENDDO ! I
00150 !
00151       ENDIF
00152 !
00153 !-----------------------------------------------------------------------
00154 !
00155       RETURN
00156       END

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