cell_area.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	cell_area

Function/Subroutine Documentation

cell_area()

subroutine cell_area

(

)

Definition at line 52 of file cell_area.f90.

!==============================================================================!
   USE all_vars
   USE mod_utils
   USE mod_par
   USE mod_spherical
   IMPLICIT NONE
   REAL(SP), ALLOCATABLE :: XX(:),YY(:) 
   REAL(SP) :: ARTMAX,ARTTOT,ARTMIN
   REAL(SP) :: ART1MAX,ART1TOT,ART1MIN
   INTEGER  :: I,J,II,J1,J2,MAX_NBRE
   
   REAL(SP) :: SBUF
   INTEGER  :: IERR
   CHARACTER(LEN=10) :: TSTR
   CHARACTER(LEN=80) ::MSG

!==============================================================================!

    IF (dbg_set(dbg_sbr)) WRITE(ipt,*) "START: CELL_AREA"
    
!---------------INITIALIZE ARRAYS----------------------------------------------!

   art  = 0.0_sp ; art1 = 0.0_sp ; art2 = 0.0_sp
   max_nbre = maxval(ntve)+1
   ALLOCATE(xx(2*max_nbre+1),yy(2*max_nbre+1))
   xx = 0.0_sp ; yy = 0.0_sp

!---------------COMPUTE AREA OF TRIANGLES USING CROSS PRODUCT------------------!
  
    DO i=1,nt 
    art(i) = (vx(nv(i,2)) - vx(nv(i,1))) * (vy(nv(i,3)) - vy(nv(i,1))) - &
             (vx(nv(i,3)) - vx(nv(i,1))) * (vy(nv(i,2)) - vy(nv(i,1)))
    END DO
   art    = abs(.5_sp*art)

!---------------COMPUTE MESH STATISTICS----------------------------------------!

   artmin = minval(art(1:n))
   artmax = maxval(art(1:n))
   arttot =    sum(art(1:n))


   IF(artmin .LT. 1.0e-6_sp) THEN
      msg = ""
      IF (par) THEN
         msg = "Proc#"
         WRITE(tstr,'(I3)') myid
         msg=trim(msg)//trim(tstr)//"; "
      END IF

      msg = trim(msg)//"Min Triangle Area="
      WRITE(tstr,'(F9.6)') artmin
      msg=trim(msg)//trim(tstr)
      
      i = minloc(art(1:n),dim=1)

      msg = trim(msg)//"; EGID="
      WRITE(tstr,'(I7)') egid(i)
      msg=trim(msg)//trim(tstr)

      WRITE(ipt,*) "*****************************"
      WRITE(ipt,*) trim(msg)
      WRITE(ipt,*) "*****************************"
   END IF



   IF (dbg_set(dbg_scl)) THEN
      WRITE(ipt,*) "!     Minimum Triangle Area: ", artmin
      WRITE(ipt,*) "!     Maximum Triangle Area: ", artmax
      WRITE(ipt,*) "!     Total Triangle Area  : ", arttot
   END IF
   IF(artmin.LT. 1.0e-6_sp) CALL warning("CELL_AREA: TRIANGLE AREA IS SMALL (LT 1e-6)")

!-------COMPUTE CONTROL VOLUME ART1: CV FOR FLUXES OF NODAL BASED VALUES-------!

   DO i=1,m
     IF(isonb(i) == 0) THEN
       DO j=1,ntve(i)
         ii=nbve(i,j)
         j1=nbvt(i,j)
         j2=j1+1-int((j1+1)/4)*3
         xx(2*j-1)=(vx(nv(ii,j1))+vx(nv(ii,j2)))*0.5_sp-vx(i)
         yy(2*j-1)=(vy(nv(ii,j1))+vy(nv(ii,j2)))*0.5_sp-vy(i)
         xx(2*j)=xc(ii)-vx(i)
         yy(2*j)=yc(ii)-vy(i)
       END DO
       xx(2*ntve(i)+1)=xx(1)
       yy(2*ntve(i)+1)=yy(1)

       DO j=1,2*ntve(i)
          art1(i)=art1(i)+0.5_sp*(xx(j+1)*yy(j)-xx(j)*yy(j+1))
       END DO
       art1(i)=abs(art1(i))
     ELSE
       DO j=1,ntve(i)
         ii=nbve(i,j)
         j1=nbvt(i,j)
         j2=j1+1-int((j1+1)/4)*3
         xx(2*j-1)=(vx(nv(ii,j1))+vx(nv(ii,j2)))*0.5_sp-vx(i)
         yy(2*j-1)=(vy(nv(ii,j1))+vy(nv(ii,j2)))*0.5_sp-vy(i)
         xx(2*j)=xc(ii)-vx(i)
         yy(2*j)=yc(ii)-vy(i)
       END DO
       j=ntve(i)+1
       ii=nbve(i,j-1)
       j1=nbvt(i,ntve(i))
       j2=j1+2-int((j1+2)/4)*3
       xx(2*j-1)=(vx(nv(ii,j1))+vx(nv(ii,j2)))*0.5_sp-vx(i)
       yy(2*j-1)=(vy(nv(ii,j1))+vy(nv(ii,j2)))*0.5_sp-vy(i)

       xx(2*j)=vx(i)-vx(i)
       yy(2*j)=vy(i)-vy(i)

       xx(2*j+1)=xx(1)
       yy(2*j+1)=yy(1)

       DO j=1,2*ntve(i)+2
        art1(i)=art1(i)+0.5_sp*(xx(j+1)*yy(j)-xx(j)*yy(j+1))
       END DO
       art1(i)=abs(art1(i))
     END IF
   ENDDO

!---------------COMPUTE MESH STATISTICS----------------------------------------!

   art1min = minval(art1(1:m))
   art1max = maxval(art1(1:m))
   art1tot =    sum(art1(1:m))

   IF(art1min .LT. 1.0e-6_sp) THEN
      msg = ""
      IF (par) THEN
         msg = "Proc#"
         WRITE(tstr,'(I3)') myid
         msg=trim(msg)//trim(tstr)//"; "
      END IF

      msg = trim(msg)//"Min Control Volume Area="
      WRITE(tstr,'(F9.6)') art1min
      msg=trim(msg)//trim(tstr)
      
      i = minloc(art1(1:m),dim=1)

      msg = trim(msg)//"; NGID="
      WRITE(tstr,'(I7)') ngid(i)
      msg=trim(msg)//trim(tstr)

      IF(isonb(i)==0) THEN
         msg = trim(msg)//"; Node is interior"
      ELSEIF(isonb(i)==1) THEN
         msg = trim(msg)//"; Node is on solid boundary"
      ELSEIF(isonb(i)==2) THEN
         msg = trim(msg)//"; Node is on open boundary"
      ELSE
         msg = trim(msg)//"; ISONB has bad value!"
      END IF


      WRITE(ipt,*) "*****************************"
      WRITE(ipt,*) trim(msg)
      WRITE(ipt,*) "*****************************"
   END IF



   IF (dbg_set(dbg_scl)) THEN
      WRITE(ipt,*) "!     Minimum Node Control Volume Area: ", art1min
      WRITE(ipt,*) "!     Maximum Node Control Volume Area: ", art1max
      WRITE(ipt,*) "!     Total Node Control Volume Area  : ", art1tot
   END IF
   IF(art1min.LT. 1.0e-6_sp) CALL warning(" CELL_AREA: NODAL CONTROL VOLUME IS SMALL (LT 1e-6)")




!---COMPUTE AREA OF CONTROL VOLUME ART2(I) = SUM(ALL TRIS SURROUNDING NODE I)--!

   DO i=1,m
     art2(i) = sum(art(nbve(i,1:ntve(i))))
   END DO
   
   art(0) = art(1) 
   art1(0) = art1(1) 
   !   IF(NT > N)ART(N+1:NT) = ART(N)
   IF(mt > m)art2(m+1:mt) = art2(m)
   IF(mt > m)art1(m+1:mt) = art1(m)
   DEALLOCATE(xx,yy)
   
   ! NOTES: SHOULD MAKE AN ARRAY TO STORE 1/ART, 1/ART2 and 1/ART2
   ! IT is faster and safer

   IF (dbg_set(dbg_log)) WRITE(ipt,*) "!  CELL AREA             :    COMPLETE"

   IF (dbg_set(dbg_sbr)) WRITE(ipt,*) "END: CELL_AREA"
   RETURN

Here is the call graph for this function:

Here is the caller graph for this function: