gradz.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\gradz.f
00002 !
00072                      SUBROUTINE GRADZ
00073 !                    ****************
00074 !
00075      &(NS,NT,NSEG,NU,NUBO,X,Y,AIRT,AIRS,CMI,JV,
00076      & ZF,DPX,DPY,DSZ,BETA,AIRST,DXIZ,DYIZ,DSP,DSM,CORR,
00077      & ELTSEG,IFABOR,MESH)
00078 !
00079 !***********************************************************************
00080 ! TELEMAC2D   V6P3                                   21/06/2013
00081 !***********************************************************************
00082 !
00083 !
00084 !
00085 !
00086 !
00087 !
00088 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00089 !| AIRS           |-->| CELL'S AREAS
00090 !| AIRST          |-->| AREAS OF SUBTRIANGLES WITHIN THE CELLS
00091 !| AIRT           |-->| TRIANGLES' AREAS
00092 !| BETA           |-->|  EXTRAPOLATION COEFFICIENT
00093 !| CMI            |-->| COORDINATES OF THE INTERFACE MIDDLE POINT
00094 !| DPX,DPY        |-->| GRADIENT OF P1 BASE FUNCTIONS
00095 !|                |   | PER TRIANGLE
00096 !| DSM            |<->| EXTRAPOLATED GRADIENTS
00097 !| CORR           |<->| CORRECTION TO HAVE CONSERVATION
00098 !| DSZ            |<--| VARIATION OF Z FOR ORDRE 2
00099 !| DXIZ,DYIZ,DSP  |<->| WORKING TABLES
00100 !| JV             |-->| NUMBER OF TRIANGLE IN WHICH IS LOCATED
00101 !|                |   | THE INTERFACE MIDDLE POINT
00102 !| NS             |-->| TOTAL NUMER OF NODES IN THE MESH
00103 !| NSEG           |-->| TOTAL NUMER OF SEGMENTS IN THE MESH
00104 !| NT             |-->| TOTAL NUMBER OF ELEMENTS IN THE MESH
00105 !| NU             |-->| NUMBERING OF NODES IN THE TRIANGLE
00106 !| NUBO           |-->| NUMBERS OF THE TWO NODES FORMING ONE EDGE (SEGMENT)
00107 !| X,Y            |-->| COORDINATES IF THE NODES
00108 !| ZF             |-->| BATHYMETRY
00109 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00110 !
00111       USE BIEF_DEF
00112       USE INTERFACE_TELEMAC2D, EX_GRADZ => GRADZ
00113       IMPLICIT NONE
00114       INTEGER LNG,LU
00115       COMMON/INFO/LNG,LU
00116 !
00117 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00118 !
00119       INTEGER, INTENT(IN)             :: NS,NT,NSEG
00120       INTEGER, INTENT(IN)             :: NU(NT,3),NUBO(2,NSEG),JV(*)
00121       DOUBLE PRECISION, INTENT(INOUT) :: DSZ(2,*)
00122       INTEGER, INTENT(IN)             :: IFABOR(NT,3)
00123       INTEGER, INTENT(IN)             :: ELTSEG(NT,3)
00124       DOUBLE PRECISION, INTENT(IN)    :: X(NS),Y(NS),AIRT(NT),AIRS(NS)
00125       DOUBLE PRECISION, INTENT(INOUT) :: DXIZ(NS),DYIZ(NS)
00126       DOUBLE PRECISION, INTENT(INOUT) :: DSP(NS),DSM(NS),CORR(NS),BETA
00127       DOUBLE PRECISION, INTENT(IN)    :: DPX(3,NT),DPY(3,NT)
00128       DOUBLE PRECISION, INTENT(IN)    :: CMI(2,*),AIRST(2,*),ZF(NS)
00129       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH
00130 !
00131 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00132 !
00133       INTEGER IS,I1,I2,I3,JT,J,NSG,NUBO1,NUBO2,ILIM,I
00134       DOUBLE PRECISION AIRJ,DXTZ,DYTZ,AIX,AIY,AJX,AJY,FACT,TEMPOR
00135       DOUBLE PRECISION GRADI,GRADJ,GRIJ,GRJI,AMDS,DSH
00136 !
00137 !-----------------------------------------------------------------------
00138 !
00139       ! MINMOD LIMITER
00140       ! EXTRAPOLATION COEFFICIENT
00141       ILIM = 1
00142       BETA = 1.D0
00143 ! INITIALISATION
00144       CALL OV( 'X=C     ' ,DXIZ   ,DXIZ  ,DXIZ  ,0.D0,NS)
00145       CALL OV( 'X=C     ' ,DYIZ   ,DYIZ  ,DYIZ  ,0.D0,NS)
00146       CALL OV( 'X=C     ' ,DSP    ,DSP   ,DSP   ,0.D0,NS)
00147       CALL OV( 'X=C     ' ,DSM    ,DSM   ,DSM   ,0.D0,NS)
00148       DSZ(1,1:NSEG)=(/(0.D0,I=1,NSEG)/)
00149       DSZ(2,1:NSEG)=(/(0.D0,I=1,NSEG)/)
00150 !
00151 !************************************************************************
00152 !     THIS LOOP IS TO COMPUTE GRAD(Z) WITH EQUATION 5.2
00153 !************************************************************************
00154       DO JT=1,NT
00155 !
00156         I1 = NU(JT,1)
00157         I2 = NU(JT,2)
00158         I3 = NU(JT,3)
00159 !
00160         AIRJ = AIRT(JT)
00161         DXTZ = ZF(I1)*DPX(1,JT)+ZF(I2)*DPX(2,JT)+ZF(I3)*DPX(3,JT) ! GRAD_X(Z)|_Tk
00162         DYTZ = ZF(I1)*DPY(1,JT)+ZF(I2)*DPY(2,JT)+ZF(I3)*DPY(3,JT) ! GRAD_Y(Z)|_Tk
00163 !
00164         TEMPOR   = AIRJ*DXTZ
00165         DXIZ(I1) = DXIZ(I1) + TEMPOR ! SUM( |C_k|*GRAD_X(Z)|_Tk )
00166         DXIZ(I2) = DXIZ(I2) + TEMPOR ! SAME AS I1
00167         DXIZ(I3) = DXIZ(I3) + TEMPOR ! SAME AS I1 AND I2
00168 !
00169         TEMPOR   = AIRJ*DYTZ
00170         DYIZ(I1) = DYIZ(I1) + TEMPOR ! SUM( |C_k|*GRAD_Y(Z)|_Tk )
00171         DYIZ(I2) = DYIZ(I2) + TEMPOR ! SAME AS I1
00172         DYIZ(I3) = DYIZ(I3) + TEMPOR ! SAME AS I1 AND I2
00173 !
00174       ENDDO
00175 !     FOR PARALLELILSM
00176       IF(NCSIZE.GT.1)THEN          ! NPON,NPLAN,ICOM,IAN
00177         CALL PARCOM2(DXIZ,DYIZ,DYIZ,NS,1,2,2,MESH )
00178       ENDIF
00179 !
00180       DO IS=1,NS
00181         DXIZ(IS) = DXIZ(IS)/(3.D0*AIRS(IS))  ! DIVIDE BY SUM(|C_k|)
00182         DYIZ(IS) = DYIZ(IS)/(3.D0*AIRS(IS))  ! DIVIDE BY SUM(|C_k|)
00183       ENDDO
00184 ! ************************************************************************
00185 !  GRAD(Z)_i IS NOW BUILT BY EQUATION 5.2
00186 ! ************************************************************************
00187 !    REBUILDS BY INTERFACE
00188 !
00189       DO NSG=1,NSEG
00190         FACT  = AIRST(1,NSG)! USEFUL FOR PARALLELISM
00191 !
00192         J     = JV(NSG) ! THIS THE TRIANGLE IN WHICH IS LOCATED CMI
00193         IF(NCSIZE.GT.1.AND.J.EQ.0)CYCLE  ! THAT MEANS CMI IS NOT LOCATED IN TRIANGLE J
00194         IF(J.EQ.0)THEN
00195             WRITE(LU,*)'@GRADZ: PROBLEM TO RETRIEVE ELEMENT'
00196             WRITE(LU,*)'IN WHICH IS LOCATED CMI'
00197             WRITE(LU,*)'LOOKING FOR EDGE NUMBER',NSG
00198             WRITE(LU,*)'WITH NODES NUMBER',NUBO(1,NSG),NUBO(2,NSG)
00199             CALL PLANTE(1)
00200             STOP
00201         ENDIF
00202 !
00203         NUBO1 = NUBO(1,NSG)
00204         NUBO2 = NUBO(2,NSG)
00205 !
00206         AIX   = CMI(1,NSG)-X(NUBO1) ! THESE ARE COORDINATES OF
00207         AIY   = CMI(2,NSG)-Y(NUBO1) !  VECTOR PM (EQ 5.1)
00208         AJX   = CMI(1,NSG)-X(NUBO2) ! P: NUBO1 OR NUBO2
00209         AJY   = CMI(2,NSG)-Y(NUBO2) ! M: CMI(NSG)
00210 !
00211 !        NODE GRADIENTS (PM.GRAD(Z) eq 5.1 of audusse paper)
00212 !
00213         GRADI = AIX*DXIZ(NUBO1) + AIY*DYIZ(NUBO1)
00214         GRADJ = AJX*DXIZ(NUBO2) + AJY*DYIZ(NUBO2)
00215 !
00216         I1 = NU(J,1)
00217         I2 = NU(J,2)
00218         I3 = NU(J,3)
00219 !
00220 !        GRADIENT BY TRIANGLE (GRAD(W)_M=GRAD(W)|Tk)
00221 !
00222         DXTZ =ZF(I1)*DPX(1,J) +ZF(I2)*DPX(2,J) + ZF(I3)*DPX(3,J)
00223         DYTZ =ZF(I1)*DPY(1,J) +ZF(I2)*DPY(2,J) + ZF(I3)*DPY(3,J)
00224 !
00225         GRIJ  = AIX*DXTZ + AIY*DYTZ
00226         GRJI  = AJX*DXTZ + AJY*DYTZ
00227 !
00228 !    EXTRAPOLATES AND CAPS  (EQUATION 5.4)
00229 !
00230         DSZ(1,NSG)  =  EXLIM(ILIM,BETA,GRADI,GRIJ )
00231         DSZ(2,NSG)  =  EXLIM(ILIM,BETA,GRADJ,GRJI )
00232 !   FOR PARALLELILSM nomore necessary while we use icom=1 in the next parcom2
00233 !         IF(NCSIZE.GT.1)THEN
00234 !          ! SEE IF WE ARE IN HALO REGION
00235 !          IF(IFABOR(J,1).EQ.-2
00236 !     &   .OR.IFABOR(J,2).EQ.-2
00237 !     &   .OR.IFABOR(J,3).EQ.-2)THEN
00238 !            FACT=0.5D0*AIRST(1,NSG)! TO ACCOUNT ONLY FOR THE HALF OD AIRST
00239 !          ENDIF
00240 !         ENDIF
00241 !
00242         IF(DSZ(1,NSG).GE.0.D0) THEN
00243           DSP(NUBO1) = DSP(NUBO1) + FACT*DSZ(1,NSG)
00244         ELSE
00245           DSM(NUBO1) = DSM(NUBO1) - FACT*DSZ(1,NSG)
00246         ENDIF
00247         IF(DSZ(2,NSG).GE.0.) THEN
00248           DSP(NUBO2) = DSP(NUBO2) + FACT*DSZ(2,NSG)
00249         ELSE
00250           DSM(NUBO2) = DSM(NUBO2) - FACT*DSZ(2,NSG)
00251         ENDIF
00252 !
00253       ENDDO
00254 !  FOR PARALLELILSM
00255       IF(NCSIZE.GT.1)THEN      ! NPON,NPLAN,ICOM,IAN , HERE ICOM=1 VALUE WITH MAX | |
00256         CALL PARCOM2(DSP,DSM,DSM,NS,1,1,2,MESH )
00257       ENDIF
00258 !
00259 !  COMPUTES THE CORRECTIONS NECESSARY TO HAVE CONSERVATION
00260 !
00261       DO IS=1,NS
00262         CORR(IS) =  DSM(IS) - DSP(IS)
00263         AMDS =MAX(DSP(IS),DSM(IS))
00264         IF(AMDS.GT.0.D0) THEN
00265           CORR(IS) = CORR(IS)/AMDS
00266         ENDIF
00267       ENDDO
00268 !
00269       DO NSG=1,NSEG
00270 !
00271         NUBO1 = NUBO(1,NSG)
00272         NUBO2 = NUBO(2,NSG)
00273 !
00274         DSH        = DSZ(1,NSG)
00275         DSZ(1,NSG) = DSH +
00276      &               MIN(0.D0,CORR(NUBO1))*MAX(0.D0,DSH)+
00277      &               MAX(0.D0,CORR(NUBO1))*MAX(0.D0,-DSH)
00278 !
00279         DSH        = DSZ(2,NSG)
00280         DSZ(2,NSG) = DSH +
00281      &               MIN(0.D0,CORR(NUBO2))*MAX(0.D0,DSH)+
00282      &               MAX(0.D0,CORR(NUBO2))*MAX(0.D0,-DSH)
00283       ENDDO
00284 !-----------------------------------------------------------------------
00285 !
00286 !  FOR PARALLELILSM
00287       IF(NCSIZE.GT.1)THEN  !  X1,X2,X3,NSEG,NPLAN,ICOM,IAN,MESH,OPT,IELM )
00288         CALL PARCOM2_SEG(DSZ(1,1:NSEG),DSZ(2,1:NSEG),DSM,NSEG,1,2,2,
00289      &                   MESH,1,11)
00290       ENDIF
00291 
00292       RETURN
00293       END

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