      SUBROUTINE pstggjcoef(isurf,pe,pf,ph,pdi,pdr,pak,pgee, pem, &
      pq0s,pls)
!....................................................................
!
! compute the glasser et al. coef. at surface isurf
!
! aplet  10.1_r8 .97
!
 USE pstcom
 USE mtrik1
 USE newmet
 USE l22com
 USE temps
 USE r33com
 USE comggf
 IMPLICIT NONE
  integer,parameter :: r8 = SELECTED_REAL_KIND(12,100)
      INTEGER ISURF
      real*8 PE
      real*8 PF
      real*8 PH
      real*8 PDI
      real*8 PDR
      real*8 PAK
      real*8 PGEE
      real*8 PEM
      real*8 PQ0S
      real*8 PLS
      INTEGER NR
      real*8 G2
      real*8 PP
      real*8 QP
      real*8 GG
      real*8 ZNQP2
      real*8 ZVP
      real*8 ZVPP
      real*8 ZB2
      real*8 ZB2OG2
      real*8 ZB2G2
      real*8 ZOB2G2
      real*8 ZOG2
      real*8 ZOB2
      real*8 ZG2OB2
      real*8 ZB
      real*8 ZALPH
      real*8 ZLAMBD
      real*8 ZRADS
      real*8 Z

!
 INTEGER	 	 	 :: t
 REAL*8, dimension(nths)         :: zbsq
!
      nr = isurf
!
! new coefficients by aplet  7.1_r8 .93
!
!
      g2 = r*r * ga(nr)**2
      pp = ppa(nr)
      qp = qpa(nr)
      gg = ga(nr)*r
      znqp2 = (n * qpa(nr))**2
!
      zvp =  0._r8 
      zvpp =  0._r8 
      zb2 =  0._r8 
      zb2og2 =  0._r8 
      zb2g2 =  0._r8 
      zob2g2 =  0._r8 
      zog2 =  0._r8 
      zob2 =  0._r8  
      zg2ob2 =  0._r8 
      zb =  0._r8 
      zalph =  0._r8 
      zlambd =  0._r8 
      zrads =  0._r8 

zbsq(1:mth) = ( grpssq(1:mth,nr)+g2 )/ xsq(1:mth,nr)

!
      dth = twopi / mth
!


zvp    = sum( xjacob(1:mth,nr)                                )*dth

! xjprym = d xjacob / d Psi; Psi = psi/twopi

zvpp   = sum( xjprym(1:mth,nr)                                )*dth/twopi
  
zb2    = sum( xjacob(1:mth,nr)*zbsq(1:mth)                    )*dth
zb2og2 = sum( xjacob(1:mth,nr)*zbsq(1:mth)/grpssq(1:mth,nr)   )*dth
zob2g2 = sum( xjacob(1:mth,nr)/(zbsq(1:mth)*grpssq(1:mth,nr)) )*dth
zog2   = sum( xjacob(1:mth,nr)/grpssq(1:mth,nr)               )*dth
zob2   = sum( xjacob(1:mth,nr)/zbsq(1:mth)                    )*dth
zg2ob2 = sum( xjacob(1:mth,nr)*grpssq(1:mth,nr)/zbsq(1:mth)   )*dth

      zalph  = twopi * n * zvp
      zlambd = - qp/( twopi*zvp**3 )
!
      pe   = - pp*zb2og2*(zvpp - gg*qp*zvp/zb2)/(twopi*qp**2)
      pf   = pp**2 *(zb2og2*zob2 + g2*(zb2og2*zob2g2 - zog2**2)) &
             /(twopi*qp)**2
      ph   = gg*pp*(zog2 - zvp*zb2og2/zb2)/(twopi*qp)
      pdi = pe + pf + ph -  0.25_r8 
      pdr = pe + pf + ph*ph
!
!
      pem  = zb2og2*( zg2ob2 + g2*(zob2 - zvp**2/zb2) )/zvp**2
!
! pq0s and pls are the typical growth rate factors and layer width divided by
! eta^1/3
!
      z   =  zb2/zb2og2
      pq0s = abs( (twopi**2 *znqp2 *z)/(zvp**2 *pem )  &
       )**( 1._r8 / 3._r8 )
      pls = abs( (zvp**2 *pem *z**2)/(twopi**2 *znqp2)  &
          )**( 1._r8 / 6._r8 )
!
      if( abs(pp) <  1.E-6_r8  ) then
         pak =  9.E99_r8 
      else
         pak  = (twopi**2 *zvp**2 *zlambd)**2 *z/(pem* pp**2)
      end if
      if( abs( pa(nr) ) <  1.E-6_r8  ) then
         pgee =  9.E99_r8 
      else
         pgee = zb2/(gamma*pa(nr)*pem*zvp)
      end if    
!
      return
!
      end


