vgfpsi.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE vgfpsi
!                   *****************
!
     &(res,ielm,u,v,f,dt,xmul,cflmax,t1,t2,mesh,msk,maskel)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    SAME AS VGRADF BUT WITH THE PSI SCHEME AND
!+                SUB-ITERATIONS TO REACH STABILITY.
!
!history  J-M HERVOUET (LNHE)
!+        18/02/08
!+        V5P9
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CFLMAX         |<--| MAXIMUM CFL NUMBER
!| DT             |-->| TIME-STEP
!| F              |-->| FUNCTION F
!| IELM           |-->| TYPE OF ELEMENT OF RESULT
!| MASKEL         |-->| MASKING OF ELEMENTS
!|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
!| MESH           |-->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
!| RES            |<--| RESULTING VECTOR
!| T1             |<->| WORK ARRAY IN BIEF_OBJ STRUCTURE
!| T2             |<->| WORK ARRAY IN BIEF_OBJ STRUCTURE
!| U              |-->| X-COMPONENT OF ADVECTION FIELD
!| V              |-->| Y-COMPONENT OF ADVECTION FIELD
!| XMUL           |-->| MULTIPLICATION FACTOR
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_vgfpsi => vgfpsi
!
      USE declarations_special
      USE interface_parallel, ONLY : p_max
      IMPLICIT NONE
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      INTEGER, INTENT(IN)           :: IELM
      LOGICAL, INTENT(IN)           :: MSK
!
      DOUBLE PRECISION, INTENT(IN)  :: DT,XMUL
      DOUBLE PRECISION, INTENT(OUT) :: CFLMAX
!
!     STRUCTURES
!
      TYPE(bief_obj), INTENT(IN)      :: U,V,F,MASKEL
      TYPE(bief_obj), INTENT(INOUT)   :: RES,T1,T2
      TYPE(bief_mesh), INTENT(INOUT)  :: MESH
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      INTEGER IT,NIT
      DOUBLE PRECISION DDT
!
      INTRINSIC int
!
!-----------------------------------------------------------------------
!
!     NUMBER OF NECESSARY SUB-ITERATIONS
!
      CALL cflpsi(t1,u,v,dt,ielm,mesh,msk,maskel)
      CALL maxi(cflmax,it,t1%R,bief_nbpts(ielm,mesh))
      IF(ncsize.GT.1) cflmax=p_max(cflmax)
!
      nit = int(cflmax) + 1
!     WRITE(LU,*) 'VGFPSI : NIT=',NIT,' CFLMAX=',CFLMAX
!
      IF(nit.GT.100) THEN
        WRITE(lu,901) nit
901     FORMAT(1x,'VGFPSI: ',1i6,' SUB-ITERATIONS REQUIRED FOR THE'
     &   ,/,1x,   '        PSI SCHEME. DECREASE THE TIME-STEP')
        CALL plante(1)
        stop
      ENDIF
!
      ddt = dt/nit
!
!     T1 WILL TAKE THE SUCCESSIVE VALUES OF F
      CALL os( 'X=Y     ' , x=t1 , y=f )
      IF(nit.GT.1) THEN
        CALL vector(mesh%T,'=','MASBAS          ',ielm,
     &              1.d0,f,f,f,f,f,f,mesh,msk,maskel)
        IF(ncsize.GT.1) CALL parcom(mesh%T,2,mesh)
        CALL os('X=1/Y   ',x=mesh%T,y=mesh%T,
     &          iopt=2,infini=0.d0,zero=1.d-8)
      ENDIF
!
!     LOOP ON THE SUB-ITERATIONS
!
      CALL cpstvc(f,res)
!
      DO it=1,nit
!
        IF(nit.GT.1) THEN
          IF(it.EQ.1) CALL os('X=0     ',x=res)
          CALL vector(t2,'=','VGRADF       PSI',ielm,
     &                xmul,t1,t1,t1,u,v,v,mesh,msk,maskel)
          CALL os('X=X+CY  ',x=res,y=t2,c=1.d0/nit)
          IF(ncsize.GT.1) CALL parcom(t2,2,mesh)
          CALL os('X=X+CYZ ',x=t1,y=t2,z=mesh%T,c=-ddt/xmul)
        ELSE
          CALL vector(res,'=','VGRADF       PSI',ielm,
     &                xmul,t1,t1,t1,u,v,v,mesh,msk,maskel)
        ENDIF
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END