C This subroutine can be used to apply a linearly varying pressure on 
C specified elements. The pressure varies linearly w.r.to x,y, or z coordinate
C (user specified).
C The user needs to edit the  section marked by "User input".
C NOTE: THE PRESSURE VALUES SPECIFIED ARE INCREMENTAL VALUES
C As set up,  positive slope causes increase in pressure 
C with depth (reduction in coordinate value)

      subroutine forcem(press,th1,th2,nn,n)
      implicit real*8 (a-h,o-z)                                               dp
      dimension n(7)
c* * * * * *
c
c     defined non-uniformed distributed force on an element.
c
c     press        distributed load increment magnitude
c     th1          coordinate
c     th2          coordinate
c     nn           integration point number
c     n            element number,etc
c
c* * * * * *
	include '../common/creeps'
	include '../common/concom'
	include '../common/blnk'

	include '../common/prepro'
	include '../common/arrays'
	include '../common/develp'

	include '../common/space'
	include '../common/dimen'
	include '../common/array2'

	dimension nd(4),ccnode(12)

C **** User input
C specify index of the coordinate along which the pressure varies
C if x, index=1; y, index=2; z,index=3
	index = 2
C if imethod=1, the user needs to specify min pressure and slope
C if imethod=2, user needs to specify min, max pressures, base coordinate
C where min. pressure is applied and height (depth,length) of part
	imethod = 2
C 2 methods of specifying pressure distribution
	if (imethod.eq.1) then
 	  pmin = 0.
	  slope = 10.
	  xx0 = 7.4246
	else if (imethod.eq.2) then
	  pmin = 0.
	  pmax = 100.
	  xx0 = 7.4246
	  depth = 6.9393
	endif

C **** End user input

	sum1 = 0

C Identify the internal nodes of the face with applied pressure
	if(n(2).eq.21) then
	 nd(1) = lm(1)
	 nd(2) = lm(2)
	 nd(3) = lm(3)
	 nd(4) = lm(4)
	else if(n(2).eq.25) then
	 nd(1) = lm(6)
	 nd(2) = lm(5)
	 nd(3) = lm(8)
	 nd(4) = lm(7)
	else if(n(2).eq.27) then
	 nd(1) = lm(2)
         nd(2) = lm(1)
         nd(3) = lm(5)
         nd(4) = lm(6)
        else if(n(2).eq.29) then
	 nd(1) = lm(3)
         nd(2) = lm(2)
         nd(3) = lm(6)
         nd(4) = lm(7)
        else if(n(2).eq.31) then
	 nd(1) = lm(4)
         nd(2) = lm(3)
         nd(3) = lm(7)
         nd(4) = lm(8)
        else if(n(2).eq.33) then
	 nd(1) = lm(1)
         nd(2) = lm(4)
         nd(3) = lm(8)
         nd(4) = lm(5)
	endif

C Find coordinates of internal node numbers	
	do i=1,4
	 jrdpre=0
	 call vecftc(ccnode,vars(ixord),ncrdmx,ncrd,nd(i),jrdpre,2,1)
	 sum1=sum1+ccnode(index)
	enddo

C Compute the center of the face
	center = sum1*0.25

C Compute the applied pressure
	if (inc.gt.0) then
	 if (imethod.eq.1) then
	   press = pmin+slope*(xx0-center)
c	   write(56,*)n(1),center,press
	 else if(imethod.eq.2) then
	   press = pmin + (pmax-pmin)*(xx0-center)/depth
	 endif
	endif

      return
      end