text:
*dk,interp_lg
      subroutine interp_lg(imsgin,xmsgin,cmsgin,msgtype,nwds,ierr2)
C
C#######################################################################
C
C $Log: interp_lg.f,v $
C Revision 2.00  2007/11/05 19:45:59  spchu
C Import to CVS
C
CPVCS
CPVCS       Rev 1.1   30 Sep 2004 09:15:52   dcg
CPVCS    replace calls to real( with calls to dble(
CPVCS
CPVCS       Rev 1.0   20 Jul 2000 14:03:26   bap
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.1   29 Jun 2000 10:53:54   bap
CPVCS    Put in PVCS logging info
C
      implicit none
C
      include 'local_element.h'
C
      integer nwds
      integer imsgin(nwds), msgtype(nwds)
      real*8 xmsgin(nwds)
      character*(*) cmsgin(nwds)
      integer i1,i2,i3,i4,i5,i6,i7,i8,ipointi,ipointj,i1temp,
     *  it,iflag,itetcnt,itoff,jtoff,l,nef,nen,nsdtopo,
     *  nsdgeom,leni,i,nptsl,icmotp,nef1,nen1,nsdtopo1,
     *  nsdgeom1,ier,icscode,icmotype,nelements2,npoints_new,
     *  nelements_new,mpno2,mpno1,ics,length,itype,ierr,ilen,
     *  nnodes2,ilayers,nelements,nnodes,ierror,ierr2,mbndry
      real*8 xtetvol,dy1,dz1,dsfac,dx1
C
C     __________________________________________________________________
C     Master Mesh Object.
C
      pointer (ipimt1, imt1)
      pointer (ipitp1, itp1)
      integer itp1(1000000), imt1(1000000)
      pointer (ipxic1, xic1)
      pointer (ipyic1, yic1)
      pointer (ipzic1, zic1)
      real*8 xic1(10000000), yic1(10000000), zic1(10000000)
      pointer (ipitetclr1, itetclr1(10000000))
      pointer (ipitettyp1, itettyp1(10000000))
      pointer (ipitetoff1, itetoff1(10000000))
      pointer (ipjtetoff1, jtetoff1(10000000))
      pointer (ipitet1, itet1(10000000))
      pointer (ipjtet1, jtet1(10000000))
      integer itet1,jtet1,itetclr1,itettyp1,itetoff1,jtetoff1
C
C     __________________________________________________________________
C     Slave Mesh Object.
C
      pointer (ipimt2, imt2)
      pointer (ipitp2, itp2)
      integer itp2(1000000), imt2(1000000)
      pointer (ipxic2, xic2)
      pointer (ipyic2, yic2)
      pointer (ipzic2, zic2)
      real*8 xic2(10000000), yic2(10000000), zic2(10000000)
      pointer (ipitetclr2, itetclr2(10000000))
      pointer (ipitettyp2, itettyp2(10000000))
      pointer (ipitetoff2, itetoff2(10000000))
      pointer (ipjtetoff2, jtetoff2(10000000))
      pointer (ipitet2, itet2(10000000))
      pointer (ipjtet2, jtet2(10000000))
      integer itet2,jtet2,itetclr2,itettyp2,itetoff2,jtetoff2
C
      pointer (ipisetwd, isetwd)
      integer isetwd(1000000)
C
      pointer (ipmpary1, mpary1)
      pointer (ipmpary2, mpary2)
      integer mpary1(1000000), mpary2(1000000)
C
      pointer (ipipflag1, ipflag1)
      integer ipflag1(1000000)
      pointer (ipitflag1, itflag1)
      integer itflag1(1000000)
C
      character*132 logmess
      character*32 isubname, cmo, cmonew
      character*8 cglobal, cdefault
      integer ip1,ip2,ip3,ip4,ip5,ip6
C
      real*8 xicvol(100), yicvol(100), zicvol(100)
C
C#######################################################################
C
      isubname='interp'
      cglobal='global'
      cdefault='default'
C
      if(cmsgin(3).eq.'-cmo-') then
         call cmo_get_name(cmo,ierror)
      else
         cmo = cmsgin(3)
      endif
      call cmo_exist(cmo,ierror)
      if(ierror.ne.0) then
         write(logmess,'(a)')
     &        'Error in subroutine extrude: input MO does not exist'
         call writloga('default',0,logmess,0,ierror)
         ierror = 1
         go to 9999
      endif
C
      call cmo_get_info('nnodes',cmo,nnodes,ilen,itype,ierr)
      call cmo_get_info('nelements',cmo,nelements,ilen,itype,ierr)
C
      ilayers=imsgin(5)
      length=nnodes
      call mmgetblk('mpary1',isubname,ipmpary1,length,1,ics)
      call mmgetblk('mpary2',isubname,ipmpary2,length,1,ics)
      if(msgtype(6).eq.1) then
         ip1=imsgin(6)
         ip2=imsgin(7)
         ip3=imsgin(8)
         mpno1=0
         do i1=ip1,ip2,ip3
            mpno1=mpno1+1
            mpary1(mpno1)=i1
         enddo
      else
         call cmo_get_info('isetwd',cmo,
     *                     ipisetwd,leni,icmotype,ierror)
         call cmo_get_info('itp1',cmo,ipitp1,leni,icmotype,ics)
         call pntlimc(cmsgin(6),cmsgin(7),cmsgin(8),
     *                ipmpary1,mpno1,nnodes,isetwd,itp1)
      endif
      if(msgtype(9).eq.1) then
         ip4=imsgin(9)
         ip5=imsgin(10)
         ip6=imsgin(11)
         mpno2=0
         do i1=ip4,ip5,ip6
            mpno2=mpno2+1
            mpary1(mpno2)=i1
         enddo
      else
         call cmo_get_info('isetwd',cmo,
     *                     ipisetwd,leni,icmotype,ierror)
         call cmo_get_info('itp1',cmo,ipitp1,leni,icmotype,ics)
         call pntlimc(cmsgin(9),cmsgin(10),cmsgin(11),
     *                ipmpary2,mpno2,nnodes,isetwd,itp1)
      endif
C
      cmonew=cmsgin(2)
C
      call cmo_exist(cmonew,icscode)
      if(icscode.eq.0) then
         call cmo_get_info('nnodes',cmonew,nnodes2,ilen,itype,ierr)
         call cmo_get_info('nelements',cmonew,
     *                     nelements2,ilen,itype,ierr)
         call cmo_select(cmonew,icscode)
         npoints_new=nnodes2+(2+ilayers)*mpno1
         nelements_new=nelements2+(2+ilayers-1)*nelements
      else
         nnodes2=0
         nelements2=0
         call cmo_derive(cmonew,cmo,icscode)
         npoints_new=(2+ilayers)*mpno1
         nelements_new=(2+ilayers-1)*nelements
      endif
C
      call cmo_get_info('ndimensions_geom',cmo,nsdgeom1,ilen,itype,ier)
      call cmo_get_info('ndimensions_topo',cmo,nsdtopo1,ilen,itype,ier)
      call cmo_get_info('nodes_per_element',cmo,nen1,ilen,itype,ierr)
      call cmo_get_info('faces_per_element',cmo,nef1,ilen,itype,ierr)
      if(nsdtopo1.eq.3) then
         print *,"Grid topology too high: ",nsdtopo1
         goto 9999
      endif
      if(nen1.eq.nelmnen(ifelmhyb) .and.
     *   nef1.eq.nelmnef(ifelmhyb)) then
         nsdgeom=3
         if(nsdtopo1.eq.1) then
            nsdtopo=2
         else
            nsdtopo=3
         endif
         nen=nelmnen(ifelmhyb)
         nef=nelmnef(ifelmhyb)
      else
         if(nsdtopo1.eq.2) then
            if(nen1.eq.nelmnen(ifelmtri) .and.
     *         nef1.eq.nelmnef(ifelmtri)) then
               nen=nelmnen(ifelmpri)
               nef=nelmnef(ifelmpri)
            elseif(nen1.eq.nelmnen(ifelmqud) .and.
     *             nef1.eq.nelmnef(ifelmqud)) then
               nen=nelmnen(ifelmhex)
               nef=nelmnef(ifelmhex)
            endif
         elseif(nsdtopo1.eq.1) then
            nsdgeom=3
            nsdtopo=2
            if(nen1.eq.nelmnen(ifelmlin) .and.
     *         nef1.eq.nelmnef(ifelmlin)) then
               nen=nelmnen(ifelmqud)
               nef=nelmnef(ifelmqud)
            endif
         endif
      endif
C
      call cmo_set_info('nnodes',cmonew,npoints_new,1,1,ierr)
      call cmo_set_info('nelements',cmonew,nelements_new,1,1,ierr)
      call cmo_set_info('ndimensions_geom',cmonew,nsdgeom,1,1,ier)
      call cmo_set_info('ndimensions_topo',cmonew,nsdtopo,1,1,ier)
      call cmo_set_info('nodes_per_element',cmonew,nen,1,1,ierr)
      call cmo_set_info('faces_per_element',cmonew,nef,1,1,ierr)
C
      call cmo_newlen(cmonew,ierr)
C
      call cmo_get_intinfo('mbndry',cmonew,mbndry,
     *                  leni,icmotp,ierror)
      call cmo_get_info('imt1',cmonew,ipimt2,leni,icmotp,ierr)
      call cmo_get_info('itp1',cmonew,ipitp2,leni,icmotp,ierr)
      call cmo_get_info('xic',cmonew,ipxic2,leni,icmotp,ierr)
      call cmo_get_info('yic',cmonew,ipyic2,leni,icmotp,ierr)
      call cmo_get_info('zic',cmonew,ipzic2,leni,icmotp,ierr)
      call cmo_get_info('itetclr',cmonew,
     *                  ipitetclr2,leni,icmotp,ier)
      call cmo_get_info('itettyp',cmonew,
     *                  ipitettyp2,leni,icmotp,ier)
      call cmo_get_info('itetoff',cmonew,
     *                  ipitetoff2,leni,icmotp,ier)
      call cmo_get_info('jtetoff',cmonew,
     *                  ipjtetoff2,leni,icmotp,ier)
      call cmo_get_info('itet',cmonew,ipitet2,leni,icmotp,ierr)
      call cmo_get_info('jtet',cmonew,ipjtet2,leni,icmotp,ierr)
C
      call cmo_get_info('imt1',cmo,ipimt1,leni,icmotp,ierr)
      call cmo_get_info('itp1',cmo,ipitp1,leni,icmotp,ierr)
      call cmo_get_info('xic',cmo,ipxic1,leni,icmotp,ierr)
      call cmo_get_info('yic',cmo,ipyic1,leni,icmotp,ierr)
      call cmo_get_info('zic',cmo,ipzic1,leni,icmotp,ierr)
      call cmo_get_info('itetclr',cmo,
     *                  ipitetclr1,leni,icmotp,ier)
      call cmo_get_info('itettyp',cmo,
     *                  ipitettyp1,leni,icmotp,ier)
      call cmo_get_info('itetoff',cmo,
     *                  ipitetoff1,leni,icmotp,ier)
      call cmo_get_info('jtetoff',cmo,
     *                  ipjtetoff1,leni,icmotp,ier)
      call cmo_get_info('itet',cmo,ipitet1,leni,icmotp,ierr)
      call cmo_get_info('jtet',cmo,ipjtet1,leni,icmotp,ierr)
C
      nptsl=mpno1
      do i=1,mpno1
         i1=mpary1(i)
         imt2(nnodes2+i)=imt1(i1)
         itp2(nnodes2+i)=itp1(i1)
         xic2(nnodes2+i)=xic1(i1)
         yic2(nnodes2+i)=yic1(i1)
         zic2(nnodes2+i)=zic1(i1)
      enddo
      do i=1,mpno2
         i1=mpary2(i)
         imt2(nnodes2+nptsl*(2+ilayers-1)+i)=imt1(i1)
         itp2(nnodes2+nptsl*(2+ilayers-1)+i)=itp1(i1)
         xic2(nnodes2+nptsl*(2+ilayers-1)+i)=xic1(i1)
         yic2(nnodes2+nptsl*(2+ilayers-1)+i)=yic1(i1)
         zic2(nnodes2+nptsl*(2+ilayers-1)+i)=zic1(i1)
      enddo
      if(ilayers.gt.0) then
         do l=2,(2+ilayers)-1
            dsfac=dble(l-1)/dble(ilayers+1)
            do i=1,mpno1
               i1=mpary1(i)
               dx1=xic2(nnodes2+nptsl*(2+ilayers-1)+i1)-xic2(nnodes2+i1)
               dy1=yic2(nnodes2+nptsl*(2+ilayers-1)+i1)-yic2(nnodes2+i1)
               dz1=zic2(nnodes2+nptsl*(2+ilayers-1)+i1)-zic2(nnodes2+i1)
               imt2(nnodes2+nptsl*(l-1)+i1)=imt2(nnodes2+i1)
               itp2(nnodes2+nptsl*(l-1)+i1)=itp2(nnodes2+i1)
               xic2(nnodes2+nptsl*(l-1)+i1)=xic2(nnodes2+i1)+dx1*dsfac
               yic2(nnodes2+nptsl*(l-1)+i1)=yic2(nnodes2+i1)+dy1*dsfac
               zic2(nnodes2+nptsl*(l-1)+i1)=zic2(nnodes2+i1)+dz1*dsfac
            enddo
         enddo
      endif
C
      length=nnodes
      call mmgetblk('ipflag1',isubname,ipipflag1,length,1,icscode)
      do i=1,nnodes
         ipflag1(i)=0
      enddo
      length=nelements
      call mmgetblk('itflag1',isubname,ipitflag1,length,1,icscode)
      do it=1,nelements
         itflag1(it)=0
      enddo
      do i=1,mpno1
         i1=mpary1(i)
         ipflag1(i1)=i
      enddo
      do it=1,nelements
         iflag=0
         do i=1,nelmnen(itettyp1(it))
            i1=itet1(itetoff1(it)+i)
            if(ipflag1(i1).gt.0) then
               iflag=iflag+1
            endif
         enddo
         if(iflag.eq.nelmnen(itettyp1(it))) then
            itflag1(it)=1
         endif
      enddo
C
      if(nelements2.le.0) then
         itetcnt=0
         itoff=0
         jtoff=0
      else
         itetcnt=nelements2
         itoff=itetoff2(nelements2)+nelmnen(itettyp2(nelements2))
         jtoff=jtetoff2(nelements2)+nelmnef(itettyp2(nelements2))
      endif
      do l=1,2+ilayers-1
         do it=1,nelements
            if(itflag1(it).eq.1) then
               if(itettyp1(it).eq.ifelmlin) then
                  i1=nnodes2+nptsl*(l-1)+itet1(itetoff1(it)+1)
                  i2=nnodes2+nptsl*(l-1)+itet1(itetoff1(it)+2)
                  i4=nptsl+i1
                  i3=nptsl+i2
                  itetcnt=itetcnt+1
                  itetclr2(itetcnt)=itetclr1(it)
                  itettyp2(itetcnt)=ifelmqud
                  itetoff2(itetcnt)=itoff
                  jtetoff2(itetcnt)=jtoff
                  itoff=itoff+nelmnen(itettyp2(itetcnt))
                  jtoff=jtoff+nelmnef(itettyp2(itetcnt))
                  itet2(itetoff2(itetcnt)+1)=i1
                  itet2(itetoff2(itetcnt)+2)=i2
                  itet2(itetoff2(itetcnt)+3)=i3
                  itet2(itetoff2(itetcnt)+4)=i4
                  do i=1,nelmnef(itettyp2(it))
                     jtet2(jtetoff2(itetcnt)+i)=-1
                  enddo
                  do i=1,nelmnen(itettyp2(itetcnt))
                     i1temp=itet2(itetoff2(itetcnt)+i)
                     xicvol(i)=xic2(i1temp)
                     yicvol(i)=yic2(i1temp)
                     zicvol(i)=zic2(i1temp)
                  enddo
                  call volume_element(itettyp2(itetcnt),
     *                                xicvol,yicvol,zicvol,
     *                                xtetvol)
                  if(xtetvol.le.0.0d+00) then
                     itet2(itetoff2(itetcnt)+1)=i1
                     itet2(itetoff2(itetcnt)+2)=i4
                     itet2(itetoff2(itetcnt)+3)=i3
                     itet2(itetoff2(itetcnt)+4)=i2
                  endif
               elseif(itettyp1(it).eq.ifelmtri) then
                  i1=nnodes2+nptsl*(l-1)+itet1(itetoff1(it)+1)
                  i2=nnodes2+nptsl*(l-1)+itet1(itetoff1(it)+2)
                  i3=nnodes2+nptsl*(l-1)+itet1(itetoff1(it)+3)
                  i4=nptsl+i1
                  i5=nptsl+i2
                  i6=nptsl+i3
                  itetcnt=itetcnt+1
                  itetclr2(itetcnt)=itetclr1(it)
                  itettyp2(itetcnt)=ifelmpri
                  itetoff2(itetcnt)=itoff
                  jtetoff2(itetcnt)=jtoff
                  itoff=itoff+nelmnen(itettyp2(itetcnt))
                  jtoff=jtoff+nelmnef(itettyp2(itetcnt))
                  itet2(itetoff2(itetcnt)+1)=i1
                  itet2(itetoff2(itetcnt)+2)=i2
                  itet2(itetoff2(itetcnt)+3)=i3
                  itet2(itetoff2(itetcnt)+4)=i4
                  itet2(itetoff2(itetcnt)+5)=i5
                  itet2(itetoff2(itetcnt)+6)=i6
                  do i=1,nelmnef(itettyp2(it))
                     jtet2(jtetoff2(itetcnt)+i)=-1
                  enddo
                  do i=1,nelmnen(itettyp2(itetcnt))
                     i1temp=itet2(itetoff2(itetcnt)+i)
                     xicvol(i)=xic2(i1temp)
                     yicvol(i)=yic2(i1temp)
                     zicvol(i)=zic2(i1temp)
                  enddo
                  call volume_element(itettyp2(itetcnt),
     *                                xicvol,yicvol,zicvol,
     *                                xtetvol)
                  if(xtetvol.le.0.0d+00) then
C*****                     itet2(itetoff2(itetcnt)+1)=i1
C*****                     itet2(itetoff2(itetcnt)+2)=i3
C*****                     itet2(itetoff2(itetcnt)+3)=i2
C*****                     itet2(itetoff2(itetcnt)+4)=i4
C*****                     itet2(itetoff2(itetcnt)+5)=i6
C*****                     itet2(itetoff2(itetcnt)+6)=i5
                  endif
               elseif(itettyp1(it).eq.ifelmqud) then
                  i1=nnodes2+nptsl*(l-1)+itet1(itetoff1(it)+1)
                  i2=nnodes2+nptsl*(l-1)+itet1(itetoff1(it)+2)
                  i3=nnodes2+nptsl*(l-1)+itet1(itetoff1(it)+3)
                  i4=nnodes2+nptsl*(l-1)+itet1(itetoff1(it)+4)
                  i5=nptsl+i1
                  i6=nptsl+i2
                  i7=nptsl+i3
                  i8=nptsl+i4
                  itetcnt=itetcnt+1
                  itetclr2(itetcnt)=itetclr1(it)
                  itettyp2(itetcnt)=ifelmhex
                  itetoff2(itetcnt)=itoff
                  jtetoff2(itetcnt)=jtoff
                  itoff=itoff+nelmnen(itettyp2(itetcnt))
                  jtoff=jtoff+nelmnef(itettyp2(itetcnt))
                  itet2(itetoff2(itetcnt)+1)=i1
                  itet2(itetoff2(itetcnt)+2)=i2
                  itet2(itetoff2(itetcnt)+3)=i3
                  itet2(itetoff2(itetcnt)+4)=i4
                  itet2(itetoff2(itetcnt)+5)=i5
                  itet2(itetoff2(itetcnt)+6)=i6
                  itet2(itetoff2(itetcnt)+7)=i7
                  itet2(itetoff2(itetcnt)+8)=i8
                  do i=1,nelmnef(itettyp2(it))
                     jtet2(jtetoff2(itetcnt)+i)=-1
                  enddo
                  do i=1,nelmnen(itettyp2(itetcnt))
                     i1temp=itet2(itetoff2(itetcnt)+i)
                     xicvol(i)=xic2(i1temp)
                     yicvol(i)=yic2(i1temp)
                     zicvol(i)=zic2(i1temp)
                  enddo
                  call volume_element(itettyp2(itetcnt),
     *                                xicvol,yicvol,zicvol,
     *                                xtetvol)
                  if(xtetvol.le.0.0d+00) then
                     itet2(itetoff2(itetcnt)+1)=i1
                     itet2(itetoff2(itetcnt)+2)=i4
                     itet2(itetoff2(itetcnt)+3)=i3
                     itet2(itetoff2(itetcnt)+4)=i2
                     itet2(itetoff2(itetcnt)+5)=i5
                     itet2(itetoff2(itetcnt)+6)=i8
                     itet2(itetoff2(itetcnt)+7)=i7
                     itet2(itetoff2(itetcnt)+8)=i6
                     do i=1,nelmnen(itettyp2(itetcnt))
                        i1temp=itet2(itetoff2(itetcnt)+i)
                        xicvol(i)=xic2(i1temp)
                        yicvol(i)=yic2(i1temp)
                        zicvol(i)=zic2(i1temp)
                     enddo
                     call volume_element(itettyp2(itetcnt),
     *                                   xicvol,yicvol,zicvol,
     *                                   xtetvol)
                  endif
               else
                  print *,"Illegal element to extrude (interp): ",it
               endif
            endif
         enddo
      enddo
C
      call cmo_set_info('nelements',cmonew,itetcnt,1,1,ierr)
C
      ipointi=nnodes2+1
      ipointj=npoints_new
      call cmo_set_info('ipointi',cmonew,
     *                ipointi,1,1,icscode)
      if (icscode .ne. 0) call x3d_error(isubname,'cmo_get_info')
      call cmo_set_info('ipointj',cmonew,
     *                ipointj,1,1,icscode)
      if (icscode .ne. 0) call x3d_error(isubname,'cmo_get_info')
C
      call geniee_cmo(cmonew)
C
      goto 9999
 9999 continue
C
      call mmrelprt(isubname,icscode)
C
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


text:
       subroutine core
*
*
*     calculate core parameters: rc,vc,nc,smc,roc
*     -------------------------------------------
*
*
      include 'common.h'
*
      real*8 sx,v2x
*
      integer i,n,nx,im,ncut
*
*
      n = nt
*
*       determination number of stars to calculate the central parameters
*
      ncut = 0.01*n
      if(ncor.le.ncut) then
        nx = ncor
      else
        if(ncut.gt.nmin) then
          nx = ncut
        else
          nx = nmin
        endif
      endif
*
*
*       determination of the central parameters
*
      sx = 0.0d0
      v2x = 0.0d0
*
      do 10 i=1,nx
         im = iname(i)
         sx = sx + body(im)
   10    v2x = v2x + body(im)*(vr(im)*vr(im) + vt(im)*vt(im))
*
      v2x = v2x/sx
      roc = 3.0d0*sx/(4.0d0*pi*r(nx)**3)
      rc = 3.0d0*v2x/(4.0d0*pi*roc)
      vc = sqrt(v2x)
      rc = sqrt(rc)
*
      i = 0
      sx = 0.0d0
*
   20 continue
      i = i + 1
      im = iname(i)
      sx = sx + body(im)
      if(rc.gt.r(i)) go to 20
      smc = sx - body(im)
      nc = i-1
*
*
      return
*
      end
*
*
*
*

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
+incdir+${MOR1KX_SOC}/ve/mor1kx_soc_alt_uvm/tb
+incdir+${MOR1KX_SOC}/ve/mor1kx_soc_alt_uvm/tests

${QUARTUS_ROOTDIR}/../modelsim_ase/altera/verilog/src/altera_mf.v

-f ${MOR1KX_SOC}/rtl/rtl.f
-f ${MEMORY_PRIMITIVES}/rtl/fpga/altsyncram/rtl.f

${MOR1KX_SOC}/ve/mor1kx_soc_alt_uvm/tb/mor1kx_soc_alt_uvm_env_pkg.sv

${MOR1KX_SOC}/ve/mor1kx_soc_alt_uvm/tests/mor1kx_soc_alt_uvm_tests_pkg.sv

${MOR1KX_SOC}/ve/mor1kx_soc_alt_uvm/tb/mor1kx_soc_alt_uvm_tb.sv


module sym_cnst

  use precision, only: ip, sp, dp
  use test_vars, only: fmt_test_name, fmt_test_res

  use tridiag, only: tridiag_sym_cnst
  use tridiag_matmul, only: tridiag_matmul_sym_cnst

  implicit none

  real(sp) :: dev_sp
  real(dp) :: dev_dp

contains

  subroutine sym_cnst_tests()

    write(*, fmt_test_name, advance="no") "sym_cnst"

    dev_dp = sym_cnst_real_dp()
    write(*, fmt_test_res, advance="no") dev_dp

    dev_sp = sym_cnst_real_sp()
    write(*, fmt_test_res, advance="no") dev_sp

    dev_dp = sym_cnst_cmplx_dp()
    write(*, fmt_test_res, advance="no") dev_dp

    dev_sp = sym_cnst_cmplx_sp()
    write(*, fmt_test_res, advance="no") dev_sp

    write(*,*)

  end subroutine sym_cnst_tests

  real(dp) function sym_cnst_real_dp() result(val)
    integer(ip), parameter :: fp = dp

    include "./sym_cnst_src/real.src"
  end function sym_cnst_real_dp

  real(dp) function sym_cnst_cmplx_dp() result(val)
    integer(ip), parameter :: fp = dp

    include "./sym_cnst_src/cmplx.src"
  end function sym_cnst_cmplx_dp

  real(sp) function sym_cnst_real_sp() result(val)
    integer(ip), parameter :: fp = sp

    include "./sym_cnst_src/real.src"
  end function sym_cnst_real_sp

  real(sp) function sym_cnst_cmplx_sp() result(val)
    integer(ip), parameter :: fp = sp

    include "./sym_cnst_src/cmplx.src"
  end function sym_cnst_cmplx_sp

end module sym_cnst


text:
      program hesmf_std
      use mod_mean_esmf  ! HYCOM ESMF mean array interface
      use mod_za         ! HYCOM array I/O interface
      implicit none
c
c --- Form the std.dev. from a mean and a mean squared HYCOM ESMF archive file.
c
      character label*81,text*18,flnm*80
      logical   meansq
      integer   mntype,iweight
c
      integer          iexpt,yrflag
      double precision time_min,time_max,time_ave,time(3)
c
      call xcspmd  !define idm,jdm
      call zaiost
      lp=6
c
      iexpt  = 0
      ii     = idm
      jj     = jdm
      iorign = 1
      jorign = 1
c
c --- number of fields involved
c
      call blkini(nn,'nn    ')
c
c --- array allocation and initialiation
c
      call mean_alloc
c
c --- land masks.
c
      call getesmfd('regional.depth')
c
      call bigrid_esmf(depths)
c
c --- first the mean archive file
c
      read (*,'(a)') flnm
      write (lp,'(2a)') ' input mean file: ',flnm(1:len_trim(flnm))
      call flush(lp)
      call getesmf(flnm,time,iweight,mntype,iexpt,yrflag)
      if     (mntype.ne.1) then
        write(lp,'(/a/)') 'error not a mean file'
        stop
      endif
      time_min = time(1)
      time_max = time(2)
      time_ave = time(3)
c
      call mean_copy
c
c --- then the mean squared archive file
c
      read (*,'(a)') flnm
      write (lp,'(2a)') ' input mnsq file: ',flnm(1:len_trim(flnm))
      call flush(lp)
      call getesmf(flnm,time,iweight,mntype,iexpt,yrflag)
      if     (mntype.ne.2) then
        write(lp,'(/a/)') 'error not a mnsq file'
        stop
      endif
      if     (time_min.ne.time(1) .or.
     &        time_max.ne.time(2) .or.
     &        time_ave.ne.time(3)     ) then
        write(lp,'(/a/)') 'error mean and mnst files incompatible'
        write(lp,*) 'time_min = ',time_min,time(1)
        write(lp,*) 'time_max = ',time_max,time(2)
        write(lp,*) 'time_ave = ',time_ave,time(3)
        stop
      endif
c
c --- output the std. dev. archive.
c
      call mean_std
c
      read (*,'(a)') flnm
      write (lp,'(2a)') 'output file: ',flnm(1:len_trim(flnm))
      call flush(lp)
      mntype = 3
      call putesmf(flnm,time_min,time_max,time_ave,
     &             mntype,iexpt,yrflag)
      end

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
C->>> ---------------------------------------------> ems_clo_eta_grp <<<
c     Closes the eta group given by hdl_eta_grp by removing gaps
c     between arrays and reassigning the first free pointer in the
c     corresponding block.
c
      subroutine ems_clo_eta_grp(eta_grp, is)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSMMGR.INC'
      include 'SLAPCS.INC'
      include 'EMSMEM.INC'
      include 'EMSP.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      integer eta_grp(0:mx_eta_grp_rec_z-1), is(0:is_n_en_m1)
      integer p_v, p_ix, p_rec, p
      integer n_v, n_ix, n_rec
      integer mem_mgr_rt_cod

      n_v = eta_grp(eta_grp_os_n_v)
      n_ix = eta_grp(eta_grp_os_n_ix)
      n_rec = eta_grp(eta_grp_os_n_rec)
      eta_grp(eta_grp_os_mx_n_v) = n_v
      eta_grp(eta_grp_os_mx_n_ix) = n_ix
      eta_grp(eta_grp_os_mx_n_rec) = n_rec
      call ems_mem_mgr_g_p8(mem_mgr_rt_cod, is,
     &     eta_grp(eta_grp_os_hdl_v), p_v)
      if (mem_mgr_rt_cod .ne. mem_mgr_rt_cod_ok) goto 8800
      call ems_mem_mgr_g_p4(mem_mgr_rt_cod, is,
     &     eta_grp(eta_grp_os_hdl_ix), p_ix)
      if (mem_mgr_rt_cod .ne. mem_mgr_rt_cod_ok) goto 8800
      call ems_mem_mgr_g_p4(mem_mgr_rt_cod, is,
     &     eta_grp(eta_grp_os_hdl_rec), p_rec)
      if (mem_mgr_rt_cod .ne. mem_mgr_rt_cod_ok) goto 8800
c
c     Shift the eta indices to cover the gap after the values.
c     Update the handle for the indices.
c
      p = p_v + rl_wo_z*(1+n_v)
      if (p .lt. p_ix) then
         call ems_cp_i_a1(1+n_ix, is(p_ix), is(p))
         eta_grp(eta_grp_os_hdl_ix + hdl_os_blk_os) =
     &        eta_grp(eta_grp_os_hdl_ix + hdl_os_blk_os) - (p_ix-p)
         eta_grp(eta_grp_os_hdl_ix + hdl_os_p) = p
      end if
      p = p + i_wo_z*(1+n_ix)
c
c     Shift the eta records to cover the gap after the indices.
c     Update the handle for the records.
c
      if (p .lt. p_rec) then
         call ems_cp_i_a1(1+n_rec, is(p_rec), is(p))
         eta_grp(eta_grp_os_hdl_rec + hdl_os_blk_os) =
     &        eta_grp(eta_grp_os_hdl_rec + hdl_os_blk_os) - (p_rec-p)
         eta_grp(eta_grp_os_hdl_rec + hdl_os_p) = p
      end if
      p = p + i_wo_z*(1+n_rec)
c
c     Close the block.
c
      call ems_mem_mgr_clo_blk(mem_mgr_rt_cod, is,
     &     1+n_rec, i_wo_z, eta_grp(eta_grp_os_hdl_rec))
      if (mem_mgr_rt_cod .ge. mem_mgr_rt_lvl_serious) goto 8800
c 7000 continue
 7100 continue
      return
 8800 continue
      ems_msg_cod = ems_msg_lvl_serious
      goto 7100
      end

C->>> -------------------------------------------> ems_ope_u_eta_grp <<<
c     Opens an eta group with handle hdl_eta_grp for updates.
c
      subroutine ems_ope_u_eta_grp(hdl_eta_grp, ds, is)
      implicit none
      include 'EMSV.INC'
      include 'EMSMMGR.INC'
      include 'EMSPM.INC'
      include 'SLAPCS.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      include 'EMSMMGRI.INC'
      integer hdl_eta_grp(0:hdl_z_m1), is(0:is_n_en_m1)
      double precision ds(0:ds_n_en_m1)
      integer prev_mx_n_u
      integer n_wo, n_fr_wo
      integer cu_is_n_en, cu_rq_is_n_en, mx_rq_is_n_en
      logical no_po
c      integer rt_cod
      integer mem_mgr_rt_cod

 100  continue
      call ems_ope_eta_grp(
     &     no_po, n_wo, n_fr_wo,
     &     hdl_eta_grp,
     &     rsmi_eta_grp_ty, sto_pk_eta_v, sto_u_eta_se_ty,
     &     mx_n_u, -1, ds, is)
      if (ems_msg_cod .ge. ems_msg_lvl_serious) goto 7000

      if (no_po) then
c
c     ?? Should be able to see if any data such as other INVERTS can be
c     deleted and the memory manager re-compressed rather than returning
c     a serious error.
c
         if (lo_sto_mode .eq. 0 .or. mx_n_u .le. 1) go to 8000
         prev_mx_n_u = mx_n_u
         mx_n_u = n_fr_wo*mx_n_u/n_wo
         if (mx_n_u .ge. prev_mx_n_u) goto 8990
         if (mx_n_u .le. 1) go to 8000
         if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9000)
     &        n_si_it, mx_n_u
         call ems_msg_wr_li(warn_msg_n)
         usr_mx_n_u = mx_n_u
         goto 100
      endif
 7000 continue
 7100 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)n_si_it
      call ems_msg_wr_li(serious_msg_n)
      call ems_mem_mgr_g_rq_is_n_en(mem_mgr_rt_cod, is,
     &     n_wo, cu_rq_is_n_en, mx_rq_is_n_en)
      if (mem_mgr_rt_cod .ge. mem_mgr_rt_lvl_serious) goto 7100
      call ems_mem_mgr_rp_rq_ws_n_en(mem_mgr_rt_cod, is, -1,
     &     cu_rq_is_n_en, mx_rq_is_n_en, rl_wo_z, 'double precision')
      if (mem_mgr_rt_cod .ge. mem_mgr_rt_lvl_serious) goto 7100
      cu_is_n_en = is(ix_l_mgr_p) + 1
      is(ix_n_xa_i_wo_rq) = max(cu_rq_is_n_en-cu_is_n_en, 1)
      goto 7100
 8990 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9899)n_wo, n_fr_wo,
     &     prev_mx_n_u, mx_n_u
      call ems_msg_wr_li(bug_msg_n)
      goto 7000
 9000 format('Iteration ', i7,
     &     ': Reducing the maximum number of updates to ', i3,
     &     ' in order to conserve space')
 9800 format('Iteration ', i7,
     &     ': Insufficient space for UPDATE')
 9899 format('Error adjusting mx_n_u: ', 4(2x, i9))
      end

C->>> ---------------------------------------------> ems_ope_eta_grp <<<
c     Returns the handle of a new eta group with space for full length
c     or packed values according to sto_eta_v
c
c     Allows space for:
c     usr_mx_n_eta etas (if >=0) or mx_n_u.
c     usr_mx_n_el  entries (if >=0) or mx_n_eta*n_r*f(av_eta_dse).
c
c     If the values are stored in long vectors then mx_n_el refers to
c     the space allocated to the indices.
c
      subroutine ems_ope_eta_grp(
     &     no_po, n_wo, n_fr_wo,
     &     hdl_eta_grp,
     &     eta_grp_ty, sto_eta_v, sto_eta_se_ty,
     &     usr_mx_n_eta, usr_mx_n_el, ds, is)
      implicit none
      include 'EMSV.INC'
      include 'EMSMMGR.INC'
      include 'EMSPM.INC'
      include 'SLAPCS.INC'
      include 'EMSMEM.INC'
      include 'EMSP.INC'
      include 'ICTVR.INC'
      include 'RLCTVR.INC'
      include 'EMSMSG.INC'
      include 'EMSMMGRI.INC'
      logical no_po
      integer n_wo, n_fr_wo
      integer hdl_eta_grp(0:hdl_z_m1)
      integer eta_grp_ty, sto_eta_v, sto_eta_se_ty
      integer usr_mx_n_eta, usr_mx_n_el
      integer is(0:is_n_en_m1)
      double precision ds(0:ds_n_en_m1)
      integer p_eta_v0, p_eta_ix0
      integer mx_n_eta, mx_n_el, blk_n, os
      integer p_eta_grp, p_eta_se_rec
      integer mem_mgr_rt_cod

      if (usr_mx_n_eta .lt. 0) then
         mx_n_eta = mx_n_u
      else
         mx_n_eta = usr_mx_n_eta
      end if
      if (usr_mx_n_el .lt. 0) then
         if (av_eta_dse .le. 0d0) goto 8010
         mx_n_el = mx_n_eta*n_r*min(one, two*av_eta_dse) + mx_n_r
      else
         mx_n_el = usr_mx_n_el
      end if
      if (sto_eta_v .eq. sto_pk_eta_v .and. (
     &     sto_eta_se_ty .eq. full_c_eta_se_ty .or.
     &     sto_eta_se_ty .eq. full_r_eta_se_ty)) then
         if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9400)
         call ems_msg_wr_li(warn_msg_n)
      endif
      call ems_g_eta_grp_n_wo(
     &     eta_grp_ty, sto_eta_v, mx_n_eta, mx_n_el, n_wo)
      if (ems_msg_cod .ge. ems_msg_lvl_serious) go to 7000
c      if (usr_mx_n_el .ge. 0) then
c         write(*, '(a, 3i9, 2x, f5.2)')
c     &        '####INVERT####: ', mx_n_eta, mx_n_el, n_wo
c      else
c         write(*, '(a, 3i9, 2x, f5.2)')
c     &        '####UPDATE####: ', mx_n_eta, mx_n_el, n_wo
c      endif
      n_fr_wo = is(ix_n_fr_wo)
      no_po = n_wo .gt. n_fr_wo
      if (no_po) goto 7000
c
c     Open a new block to store the eta group
c
      call ems_mem_mgr_ope_blk(mem_mgr_rt_cod, is,
     &     n_wo, ope_blk_at_f_fr_p,
     &     cu_ml_n, eta_grp_blk_id, blk_n)
      if (mem_mgr_rt_cod .ne. mem_mgr_rt_cod_ok) then
         if (mem_mgr_rt_cod .eq. mem_mgr_rt_cod_serious_no_po) then
            ems_msg_cod = ems_msg_lvl_serious
            go to 7000
         else if (mem_mgr_rt_cod .ge. mem_mgr_rt_lvl_serious) then
             ems_msg_cod = ems_msg_lvl_serious
            go to 7000
         endif
      endif
c
c     Let the handle for the eta group account for the whole block.
c
      call ems_mem_mgr_nw_hdl(mem_mgr_rt_cod, is,
     &     blk_n, n_wo, i_wo_z, hdl_eta_grp)
      if (ems_msg_cod .ge. ems_msg_lvl_serious) go to 7000
c
c     Initialise the eta group record.
c
      os = hdl_eta_grp(hdl_os_blk_os)
      call ems_mem_mgr_g_p4(mem_mgr_rt_cod, is,
     &     hdl_eta_grp, p_eta_grp)
      if (mem_mgr_rt_cod .ne. mem_mgr_rt_cod_ok) goto 8800
      call ems_iz_eta_grp_rec(
     &     eta_grp_ty,
     &     sto_eta_v,
     &     mx_n_eta, mx_n_el,
     &     blk_n,
     &     os,
     &     is(p_eta_grp))
      call ems_mem_mgr_g_p4(mem_mgr_rt_cod, is,
     &     is(p_eta_grp+eta_grp_os_hdl_rec), p_eta_se_rec)
      if (mem_mgr_rt_cod .ne. mem_mgr_rt_cod_ok) goto 8800
      is(p_eta_se_rec) = no_eta_se_ty
      is(p_eta_se_rec+eta_se_rec_os_eta_ty) = no_eta_se_ty
      is(p_eta_grp + eta_grp_os_n_rec) = 1
c
c     Assign a value to the zero'th value and index to avoid an
c     unassigned variable violation in cz_r.
c
      call ems_mem_mgr_g_p8(mem_mgr_rt_cod, is,
     &     is(p_eta_grp + eta_grp_os_hdl_v), p_eta_v0)
      if (mem_mgr_rt_cod .ne. mem_mgr_rt_cod_ok) goto 8800
      call ems_mem_mgr_g_p4(mem_mgr_rt_cod, is,
     &     is(p_eta_grp + eta_grp_os_hdl_ix), p_eta_ix0)
      if (mem_mgr_rt_cod .ne. mem_mgr_rt_cod_ok) goto 8800
      ds(p_eta_v0) = zero
      is(p_eta_ix0) = 0
 7000 continue
 7100 continue
      return
 8010 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9802)av_eta_dse
      call ems_msg_wr_li(bug_msg_n)
CM      IF (emsol_deb .EQ. 1) THEN
C?      call ems_dump
CM      ENDIF
      go to 7000
 8800 continue
      ems_msg_cod = ems_msg_lvl_serious
      goto 7100
 9400 format('Allocating space for packed eta values',
     &     ' but eta set type is for storage of full vectors')
 9802 format('Invalid running average eta density of ', g11.4)
      end

C->>> --------------------------------------------> ems_iz_inv_bs_bt <<<
      subroutine ems_iz_inv_bs_bt(st, vr_in_r)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      integer st(0:mx_n_c+n_r), vr_in_r(0:n_r)
      integer c_n, r_n, vr_n

      do 10, c_n = 1, n_c
         st(c_n) = st(c_n) - iand(st(c_n), inv_bs_bt)
 10   continue
      do 15, r_n = mx_n_c+1, mx_n_c+n_r
         st(r_n) = st(r_n) - iand(st(r_n), inv_bs_bt)
 15   continue
      do 20, r_n = 1, n_r
         vr_n = vr_in_r(r_n)
         st(vr_n) = st(vr_n) + inv_bs_bt
 20   continue
      return
      end

C->>> ------------------------------------------> ems_iz_eta_grp_rec <<<
      subroutine ems_iz_eta_grp_rec(
     &     eta_grp_ty,
     &     sto_eta_v,
     &     mx_n_eta,
     &     mx_n_el,
     &     blk_n,
     &     os,
     &     eta_grp)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSMMGR.INC'
      include 'SLAPCS.INC'
      include 'ICTVR.INC'
      integer eta_grp_ty, sto_eta_v, mx_n_eta, mx_n_el, blk_n, os
      integer eta_grp(0:mx_eta_grp_rec_z-1)
      integer ems_eta_grp_rec_z
      integer mx_n_v, mx_n_ix, mx_n_rec

      call ems_g_eta_grp_a_z(
     &     eta_grp_ty, sto_eta_v, mx_n_eta, mx_n_el,
     &     mx_n_v, mx_n_ix, mx_n_rec)
c
c     Initialise the eta type, the maximum and actual numbers of etas
c     and entries stored.
c
      eta_grp(eta_grp_os_eta_grp_ty) = eta_grp_ty
      eta_grp(eta_grp_os_sto_eta_v) = sto_eta_v
      eta_grp(eta_grp_os_mx_n_eta) = mx_n_eta
      eta_grp(eta_grp_os_mx_n_el) = mx_n_el
      eta_grp(eta_grp_os_f_si_it_n) = 0
      eta_grp(eta_grp_os_n_eta) = 0
      eta_grp(eta_grp_os_n_v) = 0
      eta_grp(eta_grp_os_mx_n_v) = mx_n_v
      eta_grp(eta_grp_os_n_ix) = 0
      eta_grp(eta_grp_os_mx_n_ix) = mx_n_ix
      eta_grp(eta_grp_os_n_rec) = 0
      eta_grp(eta_grp_os_mx_n_rec) = mx_n_rec
      os = os + i_wo_z*ems_eta_grp_rec_z(eta_grp_ty)
c
c     Set the handle for the values.
c
      eta_grp(eta_grp_os_hdl_v + hdl_os_blk_n) = blk_n
      eta_grp(eta_grp_os_hdl_v + hdl_os_blk_os) = os
      eta_grp(eta_grp_os_hdl_v + hdl_os_wo_z) = rl_wo_z
      os = os + rl_wo_z*(1+mx_n_v+1)

c      if (eta_grp_ty .eq. parlp_eta_grp_ty) then
cc
cc     Set the handle for the real simplex iteration data.
cc
c         eta_grp(eta_grp_os_hdl_rl_si_da + hdl_os_blk_n) = blk_n
c         eta_grp(eta_grp_os_hdl_rl_si_da + hdl_os_blk_os) = os
c         eta_grp(eta_grp_os_hdl_rl_si_da + hdl_os_wo_z) = rl_wo_z
c         os = os + rl_wo_z*(1+rl_si_da_rec_z*mx_n_eta)
c      endif
c
c     Set the handle for the indices.
c
      eta_grp(eta_grp_os_hdl_ix + hdl_os_blk_n) = blk_n
      eta_grp(eta_grp_os_hdl_ix + hdl_os_blk_os) = os
      eta_grp(eta_grp_os_hdl_ix + hdl_os_wo_z) = i_wo_z
      os = os + i_wo_z*(1+mx_n_ix+1)
c
c     Set the handle for the eta_rec.
c
      eta_grp(eta_grp_os_hdl_rec + hdl_os_blk_n) = blk_n
      eta_grp(eta_grp_os_hdl_rec + hdl_os_blk_os) = os
      eta_grp(eta_grp_os_hdl_rec + hdl_os_wo_z) = i_wo_z
      os = os + i_wo_z*(1+mx_n_rec+1)
c      if (eta_grp_ty .eq. parlp_eta_grp_ty) then
cc
cc     Set the handle for the feasibility candidate row numbers integer
cc     simplex iteration data.
cc
c         eta_grp(eta_grp_os_hdl_cdd_r_ix + hdl_os_blk_n) = blk_n
c         eta_grp(eta_grp_os_hdl_cdd_r_ix + hdl_os_blk_os) = os
c         eta_grp(eta_grp_os_hdl_cdd_r_ix + hdl_os_wo_z) = i_wo_z
c         os = os + i_wo_z*(mx_n_eta*(1+n_r/10))
c         eta_grp(eta_grp_os_hdl_i_si_da + hdl_os_blk_n) = blk_n
c         eta_grp(eta_grp_os_hdl_i_si_da + hdl_os_blk_os) = os
c         eta_grp(eta_grp_os_hdl_i_si_da + hdl_os_wo_z) = i_wo_z
c         os = os + i_wo_z*(1+i_si_da_rec_z*mx_n_eta)
c      endif
      return
      end

C->>> -------------------------------------------> ems_iz_eta_se_rec <<<
c     Initialise the record for an eta set.
c
c     |    ty    |   n_eta  | sa_(n_eta+1) |  n_eta  |  ty  |
c           1          2            3      |     4   |   5  |
c     |<---eta_se_rec_z---->|<------1----->|<-eta_se_rec_z->|
c
      subroutine ems_iz_eta_se_rec(
     &     eta_se_ty,
     &     eta_grp_n_ix,
     &     eta_se_rec)
      implicit none
      include 'EMSV.INC'
      include 'EMSMMGR.INC'
      include 'EMSPM.INC'
      integer eta_se_ty, eta_grp_n_ix
      integer eta_se_rec(1:eta_se_rec_z+1+eta_se_rec_z+1)
      integer p
c
c     Initialise the type and number of etas at the start of the record.
c
      eta_se_rec(eta_se_rec_os_eta_ty) = eta_se_ty
      eta_se_rec(eta_se_rec_os_n_eta) = 0
c
c     Backward offsets are added to the pointer to the last entry in the
c     record.
c
      p = eta_se_rec_z + 1 + eta_se_rec_z
c
c     Record the start for the 1st eta.
c     Initialise the type and number of etas at the end of the record.
c
      eta_se_rec(p-eta_se_rec_z) = eta_grp_n_ix + 1
      eta_se_rec(p+eta_se_rec_bwd_os_n_eta) = 0
      eta_se_rec(p+eta_se_rec_bwd_os_eta_ty) = eta_se_ty
c
c     Initialise the first record in the next set to correspond to a
c     null set.
c
      eta_se_rec(p+eta_se_rec_os_eta_ty) = no_eta_se_ty
      return
      end

C->>> --------------------------------------------> ems_u_eta_se_rec <<<
c     Update the record for an eta set, corresponding to the addition
c     of a single eta in an eta group with eta_grp_n_ix indices.
c
c     |  ty  | n_eta   | sa_1 ... sa_(n_eta+1) |  n_eta  |  ty  | ty' |
c                                       -2     |    -1   |   0  |  1  |
c     |<-eta_se_rec_z->|<---eta_se_n_eta+1---->|<-eta_se_rec_z->|
c
c     |  ty  | n_eta+1 | sa_1 ... sa_(n_eta+2) | n_eta+1 |  ty  | ty' |
c                                       -1     |    0    |   1  |  2  |
c     |<-eta_se_rec_z->|<---eta_se_n_eta+2---->|<-eta_se_rec_z->|
c
      subroutine ems_u_eta_se_rec(
     &     eta_grp_n_ix,
     &     eta_se_rec)
      implicit none
      include 'EMSV.INC'
      include 'EMSMMGR.INC'
      include 'EMSPM.INC'
      integer eta_grp_n_ix
      integer eta_se_rec(*)
      integer eta_se_ty, eta_se_n_eta
      integer p
c      integer ix
c
c     Get the type and number of etas from the backward record.
c
      eta_se_ty = eta_se_rec(eta_se_rec_bwd_os_eta_ty)
      eta_se_n_eta = eta_se_rec(eta_se_rec_bwd_os_n_eta)
c
c     Forward offsets are added to the pointer to the last entry in the
c     previous set.
c
      p = -(eta_se_rec_z+eta_se_n_eta+1+eta_se_rec_z)
c      write(*,9000)(ix, ix=p, 1)
c      write(*,9000)(eta_se_rec(ix), ix=p, 1)
c
c     Increase the number of etas for the forward record.
c
      eta_se_rec(p+eta_se_rec_os_n_eta) = eta_se_n_eta + 1
c
c     Initialise the first record in the next set to correspond to a
c     null set.
c
      eta_se_rec(1+eta_se_rec_os_eta_ty) = no_eta_se_ty
c
c     Shift the type and the (updated) number of etas at the end of the
c     record.
c
      eta_se_rec(1+eta_se_rec_bwd_os_eta_ty) = eta_se_ty
      eta_se_rec(1+eta_se_rec_bwd_os_n_eta) = eta_se_n_eta + 1
c
c     Record the start for the (n_eta+1)'st eta.
c
      eta_se_rec(1-eta_se_rec_z) = eta_grp_n_ix + 1
c      write(*,9000)(ix, ix=p, 2)
c      write(*,9000)(eta_se_rec(ix), ix=p, 2)
c 9000 format(20(i4))

      return
      end

C->>> -------------------------------------------> ems_eta_grp_rec_z <<<
c     Returns the number of words of data in an eta group of given type.
c
      integer function ems_eta_grp_rec_z(eta_grp_ty)
      implicit none
      include 'EMSV.INC'
      include 'EMSMMGR.INC'
      include 'SLAPCS.INC'
      integer eta_grp_ty
c      if (eta_grp_ty .eq. rsmi_eta_grp_ty) then
         ems_eta_grp_rec_z = rsmi_eta_grp_rec_z
c      else if (eta_grp_ty .eq. parlp_eta_grp_ty) then
c         ems_eta_grp_rec_z = parlp_eta_grp_rec_z
c      end if
      return
      end

C->>> ------------------------------------------> ems_g_eta_grp_n_wo <<<
c     Gets the number of integer words which make up the eta group
c     record plus the eta group itself.
c
      subroutine ems_g_eta_grp_n_wo(
     &     eta_grp_ty, sto_eta_v, mx_n_eta, mx_n_el, n_wo)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSMMGR.INC'
      include 'SLAPCS.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      integer eta_grp_ty, sto_eta_v, mx_n_eta, mx_n_el, n_wo
      integer ems_eta_grp_rec_z
      integer eta_grp_rec_z
      integer mx_n_v, mx_n_ix, mx_n_rec

      call ems_g_eta_grp_a_z(
     &     eta_grp_ty, sto_eta_v, mx_n_eta, mx_n_el,
     &     mx_n_v, mx_n_ix, mx_n_rec)
      eta_grp_rec_z = ems_eta_grp_rec_z(eta_grp_ty)
      if (eta_grp_ty .eq. rsmi_eta_grp_ty) then
         n_wo = i_wo_z*(eta_grp_rec_z) +
     &         rl_wo_z*(1+mx_n_v  +1) +
     &          i_wo_z*(1+mx_n_ix +1) +
     &          i_wo_z*(1+mx_n_rec+1)
c      else if (eta_grp_ty .eq. parlp_eta_grp_ty) then
c         n_wo = i_wo_z*(eta_grp_rec_z) +
c     &         rl_wo_z*(1+mx_n_v  +1) +
c     &         rl_wo_z*(1+rl_si_da_rec_z*mx_n_eta) +
c     &          i_wo_z*(1+mx_n_ix +1) +
c     &          i_wo_z*(1+mx_n_rec+1) +
c     &          i_wo_z*(mx_n_eta*(1+n_r/10)) +
c     &          i_wo_z*(1+i_si_da_rec_z*mx_n_eta)
      else
         go to 8010
      end if
 7000 continue
      return
 8010 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9801)
     &     eta_grp_ty, 1, n_eta_grp_ty
      call ems_msg_wr_li(bug_msg_n)
CM      IF (emsol_deb .EQ. 1) THEN
C?      call ems_dump
CM      ENDIF
      go to 7000
 9801 format('Eta group type is ', i9,
     &     '. Valid range is ', i1, ' to ', i2)
      end

C->>> -------------------------------------------> ems_g_eta_grp_a_z <<<
      subroutine ems_g_eta_grp_a_z(
     &     eta_grp_ty, sto_eta_v, mx_n_eta, mx_n_el,
     &     mx_n_v, mx_n_ix, mx_n_rec)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSMMGR.INC'
      include 'ICTVR.INC'
      integer eta_grp_ty, sto_eta_v, mx_n_eta, mx_n_el
      integer mx_n_v, mx_n_ix, mx_n_rec

      if (sto_eta_v .eq. sto_full_eta_v) then
         mx_n_v = mx_n_eta*(1+n_r)
      else
         mx_n_v = mx_n_el
      endif
      mx_n_ix = mx_n_el
c
c     Make sure that there is enough space for:
c     a start for all the possible etas plus
c     two records and a start for the n+1'st eta in all possible sets
c     plus the null type record for the n+1'st set.
c
      mx_n_rec = mx_n_eta + (2*eta_se_rec_z+1)*mx_n_eta_se + 1
      return
      end

C->>> ------------------------------------------------> ems_g_r_o_vr <<<
c     Gets the r_o_vr corresponding to the vr_in_r supplied.
c
      subroutine ems_g_r_o_vr(
     &     vr_in_r,
     &     r_o_vr)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      integer vr_in_r(0:n_r), r_o_vr(0:mx_n_c+n_r)
      integer r_n

      do 10, r_n = 1, n_r
         r_o_vr(vr_in_r(r_n)) = r_n
 10   continue
      return
      end

C->>> ------------------------------------------> ems_g_vr_in_r_perm <<<
c     Gets the permutation from fm_vr_in_r to the vr_in_r of t_r_o_vr.
c
      subroutine ems_g_vr_in_r_perm(
     &     fm_vr_in_r,
     &     st,
     &     vr_in_r_perm)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      integer fm_vr_in_r(0:n_r), st(0:mx_n_c+n_r)
      integer vr_in_r_perm(0:n_r)
      integer r_n, vr_n, og_r_n
      double precision v, v1

      do 10, r_n = 1, n_r
         vr_n = fm_vr_in_r(r_n)
         og_r_n = iand(st(vr_n), mx_mx_ml_a_dim)
         vr_in_r_perm(r_n) = og_r_n
         if (iand(st(vr_n), bc_bt) .eq. 0) then
            if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9900)vr_n
            call ems_msg_wr_li(bug_msg_n)
CM      IF (emsol_deb .EQ. 1) THEN
C?            call ems_dump
CM      ENDIF
         end if
         if (og_r_n .le. 0) then
            if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9910)
     &           vr_n, st(vr_n), f_il_ix_n, og_r_n
            call ems_msg_wr_li(bug_msg_n)
CM      IF (emsol_deb .EQ. 1) THEN
C?            call ems_dump
CM      ENDIF
            v = 0d0
            v1 = 1d0/v
         end if
 10   continue
      return
 9900 format('In g_vr_in_r_perm, variable ', i9, ' is not basic')
 9910 format('In g_vr_in_r_perm, og_r_n = ', 4(2x, i9))
      end

C->>> ------------------------------------------------> ems_perm_inv <<<
c     Apply the permutation perm to the indices of the entries in the
c     eta file.
c
      subroutine ems_perm_inv(perm, is)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSMMGR.INC'
      include 'SLAPCS.INC'
      include 'EMSMEM.INC'
      include 'EMSP.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      include 'EMSMSGN.INC'
      integer perm(0:n_r), is(0:is_n_en_m1)
      integer p_eta_grp, p_ix
      integer n_el, el_n, eta_grp_n
      integer rt_cod

      do 20, eta_grp_n = 1, eta_fi_n_grp
         p_eta_grp = eta_fi_p_a((eta_grp_n-1)*eta_fi_p_a_rec_z+
     &        eta_fi_p_a_os_p_eta_grp)
         p_ix = eta_fi_p_a((eta_grp_n-1)*eta_fi_p_a_rec_z+
     &        eta_fi_p_a_os_p_eta_ix)
         if (lo_sto_mode .ne. 0) call ems_clo_eta_grp(is(p_eta_grp), is)
CM      IF (emsol_da .EQ. 1) THEN
C?c         if (n_r .gt. 40000)
C?c     &        call ems_wr_eta_grp_n_wo(is(p_eta_grp), is)
CM      ENDIF
         n_el = is(p_eta_grp+eta_grp_os_n_ix)
         do 10, el_n = 1, n_el
            is(p_ix+el_n) = perm(is(p_ix+el_n))
 10      continue
 20   continue
      if (lo_sto_mode .ne. 0) then
         call ems_g_inv_p(rt_cod, is)
         if (rt_cod .ne. 0) goto 8020
      endif
 7000 continue
      return
 8020 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9802)
      call ems_msg_wr_li(bug_msg_n)
      goto 7000
 9802 format('Error in ems_g_inv_p')
      end

max_stars_repo_name: ERGO-Code/EMSOL

max_issues_repo_name: ERGO-Code/EMSOL

max_forks_repo_name: ERGO-Code/EMSOL


text:
In the Zone Cyber Cafe is an internet cafe with 8 highspeed quadcore Powerhouse PCs, an Xbox, a snack bar, and a few hightech items for sale. The computers can be used for internet browsing, productivity programs, and PC gaming. It costs $6/hour to use a computer. The network and internet connections are fast, so latency is less of an issue for gamers. The facility is also available to host parties, LAN, and other events.

As of March 2008, they offer this deal for new all customers:  Use their highspeed computers for 30 minutes FREE. Valid photoID required.

Question: Is the promotion for Videogamers or all gamers (including the roll a d12 kind)?

20080318 18:52:37 nbsp Is the placed closed for UCD finals week? I tried going there today (Tuesday 3/18) but it was closed... Users/AnthonyLam

20080322 17:51:46 nbsp I dropped by today and looked inside  all the computers are gone. Even the drinks are gone. LCD monitors, keyboards, etc. are all removed. Sebastion if youre reading this  are you shutting this place down?

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine invcon (kop,fmc,deltd,gain)
      dimension cone(22,13), dfreq(22), ddeg(13)
      data dfreq/2.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
     1         13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 20.0, 22.0, 24.0,
     2         26.0, 28.0, 30.0/

      data ddeg/0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0,
     2          45.0, 50.0, 55.0, 90.0/
c***********************************************************************
c          I had to break apart the data statement for MS FORTRAN to compile
      data (cone(j, 1),j=1,22)/
     +-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,
     +-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0,-20.0/
      data (cone(j, 2),j=1,22)/
     + -9.0, -7.0, -8.2, -8.6, -8.8, -8.3, -7.7, -7.7, -8.3, -8.3, -8.0,
     + -8.0, -8.8,-10.0,-10.0,-10.0,-10.0,-10.0,-10.0,-10.0,-10.0,-10.0/
      data (cone(j, 3),j=1,22)/
     + -5.0, -3.0, -3.5, -3.8, -3.9, -4.0, -4.1, -3.9, -4.7, -5.1, -4.9,
     + -4.7, -6.4, -7.0, -6.5, -6.0, -6.4, -7.6, -8.9,-10.0,-10.0,-10.0/
      data (cone(j, 4),j=1,22)/
     + -2.1,  -.1,  -.7, -1.0, -1.3, -1.4, -1.5, -1.7, -2.0, -2.6, -2.6,
     + -2.9, -4.7, -4.7, -4.0, -3.9, -4.7, -5.6, -6.7, -8.0,-10.0,-10.0/
      data (cone(j, 5),j=1,22)/
     +  -.7,  1.3,   .8,   .4,   .1,  -.7, -2.2, -2.3, -1.7, -2.4, -5.1,
     + -6.7, -5.0, -3.8, -3.5, -3.7, -4.3, -5.1, -5.7, -6.7, -8.3,-10.0/
      data (cone(j, 6),j=1,22)/
     +  -.5,  1.5,  1.4,  1.2,   .3, -1.0, -4.2, -4.7, -3.1, -2.9, -6.6,
     + -6.9, -4.0, -3.2, -3.7, -4.2, -4.8, -5.3, -5.6, -6.1, -7.6, -9.0/
      data (cone(j, 7),j=1,22)/
     +  -.3,  1.7,  1.6,  1.5,  1.2, -1.6, -5.5, -7.2, -6.0, -5.0, -8.0,
     + -8.5, -3.5, -2.7, -4.0, -4.7, -5.3, -5.7, -6.0, -6.9, -8.1,-10.0/
      data (cone(j, 8),j=1,22)/
     +  -.4,  1.5,  1.5,  1.5,  1.5, -2.2, -7.0,-11.5, -7.5, -5.8,-10.4,
     +-12.0, -2.5, -2.0, -4.3, -5.3, -6.3, -6.9, -7.4, -8.0, -9.0,-10.0/
      data (cone(j, 9),j=1,22)/
     +  -.5,  1.5,  1.5,  1.5,  1.3, -3.0, -8.4,-10.6, -6.7, -5.2, -9.8,
     +-14.0, -2.4, -2.4, -5.0, -6.2, -7.0, -7.8, -8.3,-10.0,-10.0,-10.0/
      data (cone(j,10),j=1,22)/
     +  -.8,  1.2,  1.3,  1.2,   .3, -4.1,-10.5,-11.3, -5.9, -5.0, -8.5,
     +-17.0, -7.5, -6.0, -6.4, -6.9, -7.6, -8.3,-10.0,-10.0,-10.0,-10.0/
      data (cone(j,11),j=1,22)/
     + -1.5,   .5,   .4,   .1, -1.0, -5.4,-12.1,-10.3, -6.1, -5.2, -8.5,
     +-17.0,-15.0,-10.3, -8.5, -8.0, -8.1, -8.8,-10.0,-10.0,-10.0,-10.0/
      data (cone(j,12),j=1,22)/
     + -2.2,  -.2,  -.5, -1.0, -2.2, -7.3,-12.0, -8.5, -5.6, -5.2, -7.5,
     +-15.5,-20.0,-15.0,-10.5, -8.9, -8.9,-10.0,-10.0,-10.0,-10.0,-10.0/
      data (cone(j,13),j=1,22)/
     +-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,
     +-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0,-40.0/
c***********************************************************************
      xint(a,b,c)=a*(1.0-c)+b*c
      wtf(a,b,c)=(a-b)/(c-b)
      if ( kop .ne. 30) stop ' check kop '
c     kop=30,  inverted cone antenna
      freq=fmc
      delta=deltd

      do 10 iang = 2,13
      if ( delta .gt. ddeg(iang) ) go to 10
      ian = iang
      go to 20
   10 continue

   20 do 30 ifrq = 2,22
      if ( freq .gt. dfreq(ifrq)) go to 30
      ifr = ifrq
      go to 40
   30 continue

   40 fwt=wtf(freq,dfreq(ifr),dfreq(ifr-1))
      awt=wtf(delta,ddeg(ian),ddeg(ian-1))
      g1=xint(cone(ifr,ian),cone(ifr,ian-1),awt)
      g2=xint(cone(ifr-1,ian),cone(ifr-1,ian-1),awt)
      gain=xint(g1,g2,fwt)
      return
      end

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


text:
      SUBROUTINE MSFERM(PAR,IFAIL,FLAG)

*******************************************************************
* Subroutine to compute the squark pole masses/mixing angles and
* slepton masses/mixing angles (with thanks to S. Kraml).
*
* The (large tan(beta)) SUSY correction DELMB is computed here.
*
* Negative squark masses squared give IFAIL = 8.
*
* The running masses of the 3rd generation of squarks are first
* evaluated at a scale QSTSB taking all possible thresholds
* (i.e. possibly large logs) into account. Pole masses are computed
* including 1 loop (S)QCD rad. corrs.
* No electroweak rad. corrs. are included at present.
*
* Concerning the 1st and 2nd generation squarks, it is assumed that
* the input parameters PAR(15)=MQ, PAR(16)=MU, PAR(17)=MD are the
* running squark masses at the SUSY scale Q2 ~ MQ^2 ~ MU^2 ~ MD^2.
*
* The conventions for the cosines/sines CSF/SSF (F=T, B, L=Stau) of the
* 3rd generation sfermion mixing angle are
* F1 = CSF*FL + SSF*FR, F2 = CSF*FR - SSF*FL
* where FR, FL are the left/right handed weak eigenstates, and
* F1, F2 are the physical states ordered in mass (M_F1 < M_F2)
*******************************************************************

      IMPLICIT NONE

      INTEGER IFAIL,OMGFLAG,MAFLAG,MOFLAG,FLAG,PFLAG

      DOUBLE PRECISION PAR(*)
      DOUBLE PRECISION cor,Q2,QSTSB,pi
      DOUBLE PRECISION Wt,At,mstL,mstR
      DOUBLE PRECISION Wb,Ab,msbL,msbR
      DOUBLE PRECISION Wl,Atau,AMUON,Xl,mslL,mslR
      DOUBLE PRECISION ALSMZ,ALEMMZ,GF,g1,g2,S2TW
      DOUBLE PRECISION MS,MC,MB,MBP,MT,MTAU,MMUON,MZ,MW,X
      DOUBLE PRECISION MUR,MUL,MDR,MDL,MLR,MLL,MNL
      DOUBLE PRECISION MST1,MST2,MSB1,MSB2,MSL1,MSL2,MSNT
      DOUBLE PRECISION CST,CSB,CSL,MSMU1,MSMU2,MSNM,CSMU
      DOUBLE PRECISION MGL,nen,fac,DMSQUARK
      DOUBLE PRECISION RMST1,RMST2,S2T,RMSB1,RMSB2,S2B,XT,XB
      DOUBLE PRECISION M1,M2,M3,MQ3,MU3,MD3,ML3,ME3,MQ,MU,MD,ML,ME
      DOUBLE PRECISION MQ3P,MU3P,MD3P,MH1S,MH2S,MSS,B,GM
      DOUBLE PRECISION LQSTSB
      DOUBLE PRECISION RTOP,RBOT,COEF,ATP,ABP
      DOUBLE PRECISION INTEG,DELMB,C2T,C2B
      DOUBLE PRECISION G1Q,G2Q,GQ,ALSQ
      DOUBLE PRECISION ZHU,ZHD,ZS,H1Q,H2Q,TANBQ
      DOUBLE PRECISION HTQ,HBQ,MTQ,MBQ
      DOUBLE PRECISION LQ,KQ,ALQ,AKQ,MUQ,NUQ
      DOUBLE PRECISION XIFQ,XISQ,MUPQ,MSPQ,M3HQ
      DOUBLE PRECISION XI,HT2,HB2,MH1,MH2,L2,K2,HTAU2,AL,AK
      DOUBLE PRECISION g1s,g2s,g3s,HTOPS,HBOTS,HTAUS,SP,SIG1,SIG2,SIG3
      DOUBLE PRECISION F20,F21,F22,F14,F15
      DOUBLE PRECISION DPI11ST,DPI22ST,DPI11SB,DPI22SB

      COMMON/FLAGS/OMGFLAG,MAFLAG,MOFLAG
      COMMON/GAUGE/ALSMZ,ALEMMZ,GF,g1,g2,S2TW
      COMMON/SMSPEC/MS,MC,MB,MBP,MT,MTAU,MMUON,MZ,MW
      COMMON/SFSPEC/MUR,MUL,MDR,MDL,MLR,MLL,MNL,
     .      MST1,MST2,MSB1,MSB2,MSL1,MSL2,MSNT,
     .      CST,CSB,CSL,MSMU1,MSMU2,MSNM,CSMU
      COMMON/RADCOR/RMST1,RMST2,S2T,RMSB1,RMSB2,S2B,XT,XB
      COMMON/RADCOR2/MQ3P,MU3P,MD3P,ATP,ABP
      COMMON/RENSCALE/Q2
      COMMON/STSBSCALE/QSTSB
      COMMON/DELMB/DELMB
      COMMON/QGAUGE/G1Q,G2Q,GQ,ALSQ
      COMMON/QHIGGS/ZHU,ZHD,ZS,H1Q,H2Q,TANBQ
      COMMON/QQUARK/HTQ,HBQ,MTQ,MBQ
      COMMON/QPAR/LQ,KQ,ALQ,AKQ,MUQ,NUQ
      COMMON/QEXT/XIFQ,XISQ,MUPQ,MSPQ,M3HQ
      COMMON/SUSYMH/MH1S,MH2S,MSS
      COMMON/PFLAG/PFLAG
      COMMON/SUSYCOUP/g1s,g2s,g3s,HTOPS,HBOTS,HTAUS

      PI=4d0*DATAN(1d0)
      COEF=1d0/(16d0*PI**2)

!      WRITE(0,*)"CALL MSFERM"
!      WRITE(0,*)""

      AT=PAR(12)
      AB=PAR(13)
      ATAU=PAR(14)
      AMUON=PAR(25)

*   In order to run the sfermion masses from Q2 to QSTSB
*   tree level expressions for the soft
*   Higgs masses and various logarithms are needed:

      IF(MAFLAG.GE.0)THEN
* Otherwise: taken from COMMON/SUSYMH
       B=ALQ+NUQ
       GM= LQ*XIFQ+MUPQ*MUQ+M3HQ
       MH1S= -LQ**2*H2Q**2 - MUQ**2 + (MUQ*B+GM)/TANBQ
     .      + GQ/2d0*(H2Q**2-H1Q**2)
       MH2S= -LQ**2*H1Q**2 - MUQ**2 + (MUQ*B+GM)*TANBQ
     .      + GQ/2d0*(H1Q**2-H2Q**2)
      ENDIF

      M1=PAR(20)
      M2=PAR(21)
      M3=PAR(22)
      MGL=PAR(22)

      MQ3=PAR(7)
      MU3=PAR(8)
      MD3=PAR(9)
      ML3=PAR(10)
      ME3=PAR(11)
      MQ=PAR(15)
      MU=PAR(16)
      MD=PAR(17)
      ML=PAR(18)
      ME=PAR(19)

      LQSTSB=DLOG(Q2/QSTSB)

*   Running stop/sbottom masses squared and mixings:
*   On input, they are defined at the scale Q2.
*   First, they have to be evaluated at the scale QSTSB
*   (as used in the calculation of the pole masses)
*   Use 2-loop RGEs from Q2 to QSTSB in order to comply
*   with SLHA (input parameters in DR_bar without thresholds)

* CALL GETSUSYCOUP(PAR,IFAIL) in order to determine the gauge and
* Yukawa couplings in "COMMON/SUSYCOUP/g1s,g2s,g3s,HTOPS,HBOTS,HTAUS" at Q

       CALL GETSUSYCOUP(PAR,IFAIL)

       XI= g1S*(MH1S-MH2S+MQ3+2d0*MQ-2d0*(MU3+2d0*MU)+MD3+2d0*MD
     .   +ME3+2d0*ME-ML3+2d0*ML)
       HT2=HTOPS**2
       HB2=HBOTS**2
       HTAU2=HTAUS**2
       MH1=MH1S
       MH2=MH2S
       L2=PAR(1)**2
       K2=PAR(2)**2
       AL=PAR(5)
       AK=PAR(6)

      SP= HT2*(-3d0*MH1 - MQ3 + 4d0*MU3)
     .  + HB2*(3d0*MH2 - MQ3 - 2d0*MD3)
     .  + HTAU2*(MH2 + ML3 - 2d0*ME3) + L2*(MH2 - MH1)
     .  + (G1S/18d0 + 3d0/2d0*G2S + 8d0/3d0*G3S)*(MQ3+2d0*MQ)
     .  - (16d0/9d0*G1S + 16d0/3d0*G3S)*(MU3+2d0*MU)
     .  + (2d0/9d0*G1S + 8d0/3d0*G3S)*(MD3+2d0*MD)
     .  + (G1S/2d0 + 3d0/2d0*G2S)*(MH1-MH2-(ML3+2d0*ML))
     .  + 2d0*G1S*(ME3+2d0*ME)
      SIG1= G1S*(MH1 + MH2 + (MQ3+2d0*MQ)/3d0 + 8d0/3d0*(MU3+2d0*MU)
     .    + 2d0/3d0*(MD3+2d0*MD)+ ML3+2d0*ML + 2d0*(ME3+2d0*ME))
      SIG2=G2S*(MH1 + MH2 + 3d0*(MQ3+2d0*MQ) + ML3+2d0*ML)
      SIG3=G3S*(2d0*(MQ3+2d0*MQ) + MU3+2d0*MU + MD3+2d0*MD)

* Copy F20, F21, F22 from integs.f, and S -> XI,
* G1,G2,G3 -> G1S,G2S,G3S, MS -> MSS

      F20= HT2*(MH1+MQ3+MU3+AT**2) + HB2*(MH2+MQ3+MD3+AB**2)
     .       - g1S*M1**2/9d0 - 3d0*g2S*M2**2
     .       - 16d0/3d0*g3S*M3**2 + XI/6d0
     .       + COEF*(-10d0*HT2**2*(MH1+MQ3+MU3+2d0*AT**2)
     .       - 10d0*HB2**2*(MH2+MQ3+MD3+2d0*AB**2)
     .       - 10d0*HB2**4*(MH2+MQ3+MD3+2d0*AB**2)
     .       - HB2*HTAU2*(2d0*MH2+MQ3+MD3+ML3+ME3+(AB+ATAU)**2)
     .       - L2*HT2*(2d0*MH1+MH2+MSS+MQ3+MU3+(AL+AT)**2)
     .       - L2*HB2*(MH1+2d0*MH2+MSS+MQ3+MD3+(AL+AB)**2)
     .       + 4d0/3d0*G1S*HT2*(MH1+MQ3+MU3+AT**2+2d0*M1*(M1-AT))
     .       + 2d0/3d0*G1S*HB2*(MH2+MQ3+MD3+AB**2+2d0*M1*(M1-AB))
     .       + 199d0/54d0*G1S**2*M1**2 + 33d0/2d0*G2S**2*M2**2
     .     - 64d0/3d0*G3S**2*M3**2
     .     + 1d0/3d0*G1S*G2S*(M1**2+M2**2+M1*M2)
     .       + 16d0/27d0*G1S*G3S*(M1**2+M3**2+M1*M3)
     .     + 16d0*G2S*G3S*(M2**2+M3**2+M2*M3)
     .       + 1d0/3d0*G1S*SP + 1d0/18d0*G1S*SIG1
     .       + 3d0/2d0*G2S*SIG2 + 8d0/3d0*G3S*SIG3)

      F21= 2d0*HT2*(MH1+MQ3+MU3+AT**2)
     .       - 16d0/9d0*g1S*M1**2 - 16d0/3d0*g3S*M3**2 - 2d0*XI/3d0
     .       + COEF*(-16d0*HT2**2*(MH1+MQ3+MU3+2d0*AT**2)
     .       - 2d0*HT2*HB2*(MH1+MH2+2d0*MQ3+MU3+MD3+(AT+AB)**2)
     .       - 2d0*HT2*L2*(2d0*MH1+MH2+MSS+MQ3+MU3+(AT+AL)**2)
     .       - 2d0/3d0*G1S*HT2*(MH1+MQ3+MU3+AT**2+2d0*M1*(M1-AT))
     .       + 6d0*G2S*HT2*(MH1+MQ3+MU3+AT**2+2d0*M2*(M2-AT))
     .       + 1712d0/27d0*G1S**2*M1**2 - 64d0/3d0*G3S**2*M3**2
     .       + 256d0/27d0*G1S*G3S*(M1**2+M3**2+M1*M3)
     .     - 4d0/3d0*G1S*SP + 8d0/9d0*G1S*SIG1 + 8d0/3d0*G3S*SIG3)

      F22= 2d0*HB2*(MH2+MQ3+MD3+AB**2)
     .       - 4d0/9d0*g1S*M1**2 - 16d0/3d0*g3S*M3**2 + XI/3d0
     .       + COEF*(-16d0*HB2**2*(MH2+MQ3+MD3+2d0*AB**2)
     .       - 2d0*HB2*HT2*(MH1+MH2+2d0*MQ3+MU3+MD3+(AB+AT)**2)
     .       - 2d0*HB2*HTAU2*(2d0*MH2+MQ3+MD3+ML3+ME3+(AB+ATAU)**2)
     .       - 2d0*HB2*L2*(MH1+2d0*MH2+MSS+MQ3+MD3+(AB+AL)**2)
     .       + 2d0/3d0*G1S*HB2*(MH2+MQ3+MD3+AB**2+2d0*M1*(M1-AB))
     .       + 6d0*G2S*HB2*(MH2+MQ3+MD3+AB**2+2d0*M2*(M2-AB))
     .       + 404d0/27d0*G1S**2*M1**2 - 64d0/3d0*G3S**2*M3**2
     .       + 64d0/27d0*G1S*G3S*(M1**2+M3**2+M1*M3)
     .     + 2d0/3d0*G1S*SP + 2d0/9d0*G1S*SIG1 + 8d0/3d0*G3S*SIG3)

        MQ3P=MQ3-COEF*LQSTSB*F20
        MU3P=MU3-COEF*LQSTSB*F21
        MD3P=MD3-COEF*LQSTSB*F22

      F14= 6d0*HT2*AT + HB2*AB + L2*AL
     .     + 13d0/9d0*g1S*M1 + 3d0*g2S*M2 + 16d0/3d0*g3S*M3
     .     + COEF*(-44d0*HT2**2*AT - 5d0*HT2*HB2*(AT+AB)
     .     - 3d0*HT2*L2*(AT+AL) - 10d0*HB2**2*AB
     .     - HB2*HTAU2*(AB+ATAU) - 4d0*HB2*L2*(AB+AL)
     .     - 6d0*L2**2*AL - L2*HTAU2*(AL+ATAU)
     .     - 2d0*L2*K2*(AL+AK)
     .     + 2d0*g1S*HT2*(AT-M1) + 6d0*g2S*HT2*(AT-M2)
     .     + 16d0*g3S*HT2*(AT-M3) + 2d0/3d0*g1S*HB2*(AB-M1)
     .     - 2743d0/81d0*g1S**2*M1 - 5d0/3d0*g1S*g2S*(M1+M2)
     .     - 136d0/27d0*g1S*g3S*(M1+M3) - 15d0*g2S**2*M2
     .     - 8d0*g2S*g3S*(M2+M3) + 32d0/9d0*g3S**2*M3)

      F15= 6d0*HB2*AB + HT2*AT + HTAU2*ATAU + L2*AL
     .     + 7d0/9d0*g1S*M1 + 3d0*g2S*M2 + 16d0/3d0*g3S*M3
     .     + COEF*(-44d0*HB2**2*AB - 5d0*HT2*HB2*(AT+AB)
     .     - 3d0*HB2*HTAU2*(AB+ATAU) - 3d0*HB2*L2*(AB+AL)
     .     - 10d0*HT2**2*AT - 4d0*HT2*L2*(AT+AL)
     .     - 6d0*HTAU2**2*ATAU - 6d0*L2**2*AL - 2d0*L2*K2*(AL+AK)
     .     + 2d0/3d0*g1S*HB2*(AB-M1) + 6d0*g2S*HB2*(AB-M2)
     .     + 16d0*g3S*HB2*(AB-M3) + 4d0/3d0*g1S*HT2*(AT-M1)
     .     + 2d0*HTAU2*g1S*(ATAU-M1)
     .     - 1435d0/81d0*g1S**2*M1 - 5d0/3d0*g1S*g2S*(M1+M2)
     .     - 40d0/27d0*g1S*g3S*(M1+M3) - 15d0*g2S**2*M2
     .     - 8d0*g2S*g3S*(M2+M3) + 32d0/9d0*g3S**2*M3)

        ATP=AT-COEF*LQSTSB*F14
        ABP=AB-COEF*LQSTSB*F15

*   Running stop masses squared and mixings

      mstL= MQ3P + mtq**2 + (gQ/2d0-g1Q/3d0)*(h2q**2-h1q**2)
      mstR= MU3P + mtq**2 + g1Q/3d0*(h2q**2-h1q**2)
      Xt= ATP-muq/tanbq
      Wt= DSQRT( (mstL-mstR)**2 + 4d0*(Xt*mtq)**2)
      RMST1= 0.5d0*(mstL+mstR-Wt)
      RMST2= 0.5d0*(mstL+mstR+Wt)

      IF(RMST1.LE.1d1)THEN
       IFAIL=8
       nen=0d0
       IF(FLAG.EQ.0)RMST1=1d1
!      WRITE(0,*)"MST1^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ELSE
       nen= DSQRT((mstL-RMST1)**2 + (Xt*mtq)**2)
      ENDIF

      IF(RMST2.LE.1d1)THEN
       IFAIL=8
       IF(FLAG.EQ.0)RMST2=1d1
!      WRITE(0,*)"MST2^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(nen.EQ.0d0)THEN
       IF(mstL.LE.mstR)THEN
        CST= 1d0
       ELSE
        CST= 0d0
       ENDIF
      ELSE
       CST= -mtq*Xt/nen
      ENDIF
      S2T=2d0*CST*DSQRT(1d0-CST**2)

*   Running sbottom masses squared and mixings

      msbL= MQ3P + mbq**2 - (gQ/2d0-g1Q/6d0)*(h2q**2-h1q**2)
      msbR= MD3P + mbq**2 - g1Q/6d0*(h2q**2-h1q**2)
      Xb= ABP-muq*tanbq
      Wb= DSQRT((msbL-msbR)**2 + 4d0*(Xb*mbq)**2)
      RMSB1= 0.5d0*(msbL+msbR-Wb)
      RMSB2= 0.5d0*(msbL+msbR+Wb)

      IF(RMSB1.LE.1d1)THEN
       IFAIL=8
       nen=0d0
       IF(FLAG.EQ.0)RMSB1=1d1
!      WRITE(0,*)"MSB1^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ELSE
       nen= DSQRT((msbL-RMSB1)**2 + (Xb*mbq)**2)
      ENDIF

      IF(RMSB2.LE.1d1)THEN
       IFAIL=8
       IF(FLAG.EQ.0)RMSB2=1d1
!      WRITE(0,*)"MSB2^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(nen.EQ.0d0)THEN
       IF(msbL.LE.msbR)THEN
        CSB= 1d0
       ELSE
        CSB= 0d0
       ENDIF
      ELSE
       CSB= -mbq*Xb/nen
      ENDIF
      S2B=2*CSB*DSQRT(1d0-CSB**2)

*   Pole Squark Masses squared

      fac=ALSQ/(3d0*PI)
      MST1=RMST1
      MST2=RMST2
      MSB1=RMSB1
      MSB2=RMSB2

      CALL DMSQUARK_YUK(IFAIL,DPI11ST,DPI22ST,DPI11SB,DPI22SB)

      IF(MST1.GT.1d1)THEN
       MST1=RMST1-fac*DMSQUARK(1,RMST1,RMST2,S2T,mtq,MGL,QSTSB)
      MST1=MST1-DPI11ST
       MST2=RMST2-fac*DMSQUARK(2,RMST1,RMST2,S2T,mtq,MGL,QSTSB)
      MST2=MST2-DPI22ST
      ENDIF
      IF(MSB1.GT.1d1)THEN
       MSB1=RMSB1-fac*DMSQUARK(1,RMSB1,RMSB2,S2B,mbq,MGL,QSTSB)
      MSB1=MSB1-DPI11SB
       MSB2=RMSB2-fac*DMSQUARK(2,RMSB1,RMSB2,S2B,mbq,MGL,QSTSB)
      MSB2=MSB2-DPI22SB
      ENDIF

      IF(MST1.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MST1=1d2
!      WRITE(0,*)"MST1^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MST2.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MST2=1d2
!      WRITE(0,*)"MST2^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MSB1.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MSB1=1d2
!      WRITE(0,*)"MSB1^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MSB2.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MSB2=1d2
!      WRITE(0,*)"MSB2^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

*   First generation squark masses squared

      MUR=MU+g1Q/3d0*(h2q**2-h1q**2)
      MUL=MQ+(gQ/2d0-g1Q/3d0)*(h2q**2-h1q**2)
      MDR=MD-g1Q/6d0*(h2q**2-h1q**2)
      MDL=MQ+(-gQ/2d0+g1Q/6d0)*(h2q**2-h1q**2)

      IF(MUR.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MUR=1d2
!      WRITE(0,*)"MUR^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MUL.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MUL=1d2
!      WRITE(0,*)"MUL^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MDR.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MDR=1d2
!      WRITE(0,*)"MDR^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MDL.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MDL=1d2
!      WRITE(0,*)"MDL^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

*   Pole Squark Masses squared
*   On input they are defined at the scale Q2 that is assumed to not
*   too far from their actual masses
*   (Threshold effects are neglected)

      IF(MUR.GT.1d2)THEN
       X=M3**2/MUR
       IF(X.NE.1d0)THEN
        cor=1d0+3d0*X+1d0/2d0*(1d0-X)**2*DLOG((1d0-X)**2)
     .     -X**2*DLOG(X)+2d0*X*DLOG(Q2/MUR)
       ELSE
        cor=4d0+2d0*DLOG(Q2/MUR)
       ENDIF
       MUR=MUR*(1d0+2d0*fac*cor)
      ENDIF

      IF(MUL.GT.1d2)THEN
       X=M3**2/MUL
       IF(X.NE.1d0)THEN
        cor=1d0+3d0*X+1d0/2d0*(1d0-X)**2*DLOG((1d0-X)**2)
     .     -X**2*DLOG(X)+2d0*X*DLOG(Q2/MUL)
       ELSE
        cor=4d0+2d0*DLOG(Q2/MUL)
       ENDIF
       MUL=MUL*(1d0+2d0*fac*cor)
      ENDIF

      IF(MUL.GT.1d2)THEN
       X=M3**2/MDR
       IF(X.NE.1d0)THEN
        cor=1d0+3d0*X+1d0/2d0*(1d0-X)**2*DLOG((1d0-X)**2)
     .     -X**2*DLOG(X)+2d0*X*DLOG(Q2/MDR)
       ELSE
        cor=4d0+2d0*DLOG(Q2/MDR)
       ENDIF
       MDR=MDR*(1d0+2d0*fac*cor)
      ENDIF

      IF(MUL.GT.1d2)THEN
       X=M3**2/MDL
       IF(X.NE.1d0)THEN
        cor=1d0+3d0*X+1d0/2d0*(1d0-X)**2*DLOG((1d0-X)**2)
     .     -X**2*DLOG(X)+2d0*X*DLOG(Q2/MDL)
       ELSE
        cor=4d0+2d0*DLOG(Q2/MDL)
       ENDIF
       MDL=MDL*(1d0+2d0*fac*cor)
      ENDIF

      IF(MUR.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MUR=1d2
!      WRITE(0,*)"MUR^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MUL.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MUL=1d2
!      WRITE(0,*)"MUL^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MDR.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MDR=1d2
!      WRITE(0,*)"MDR^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MDL.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MDL=1d2
!      WRITE(0,*)"MDL^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MST1.GE.0d0)MST1=DSQRT(MST1)
      IF(MST2.GE.0d0)MST2=DSQRT(MST2)
      IF(MSB1.GE.0d0)MSB1=DSQRT(MSB1)
      IF(MSB2.GE.0d0)MSB2=DSQRT(MSB2)
      IF(MUL.GE.0d0)MUL=DSQRT(MUL)
      IF(MUR.GE.0d0)MUR=DSQRT(MUR)
      IF(MDL.GE.0d0)MDL=DSQRT(MDL)
      IF(MDR.GE.0d0)MDR=DSQRT(MDR)

* Calculation of the SUSY corrections to h_bot, DELMB, as in
* Carena et al., hep-ph/9912516

      C2T=1d0-S2T
      C2B=1d0-S2B

      DELMB=MUQ*TANBQ*(2d0/(3d0*PI)*ALSMZ*M3*INTEG(MSB1,MSB2,M3)
     .   +COEF*HTQ**2*AT*INTEG(MST1,MST2,MUQ)
     .   -COEF*G2Q*M2*(CST**2*INTEG(MST1,M2,MUQ)
     .   +(1d0-CST**2)*INTEG(MST2,M2,MUQ)
     .   +1d0/2d0*(CSB**2*INTEG(MSB1,M2,MUQ)
     .   +(1d0-CSB**2)*INTEG(MSB2,M2,MUQ))))

!      WRITE(0,*)"DELMB =",DELMB
!      WRITE(0,*)""

*   Stau masses squared and mixings

      mslL= ML3 + MTAU**2 - (gQ/2d0-g1Q/2d0)*(h2q**2-h1q**2)
      mslR= ME3 + MTAU**2 - g1Q/2d0*(h2q**2-h1q**2)
      Xl= Atau-MUQ*tanbq
      Wl= DSQRT((mslL-mslR)**2 + 4d0*(Xl*MTAU)**2)
      msl1= 0.5d0*(mslL+mslR-Wl)
      msl2= 0.5d0*(mslL+mslR+Wl)

      nen= DSQRT((mslL-msl1)**2 + (Xl*MTAU)**2)
      IF(nen.EQ.0d0)THEN
       IF(mslL.LE.mslR)THEN
        CSL= 1d0
       ELSE
        CSL= 0d0
       ENDIF
      ELSE
       CSL= -MTAU*Xl/nen
      ENDIF

      MSNT= ML3 + gQ/2d0*(h2q**2-h1q**2)

      IF(MSL1.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MSL1=1d2
!      WRITE(0,*)"MSL1^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MSL2.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MSL2=1d2
!      WRITE(0,*)"MSL2^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MSNT.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MSNT=1d2
!      WRITE(0,*)"MSNT^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MSL1.GE.0d0)MSL1=DSQRT(MSL1)
      IF(MSL2.GE.0d0)MSL2=DSQRT(MSL2)
      IF(MSNT.GE.0d0)MSNT=DSQRT(MSNT)

*   Smuon masses squared and mixings

      mslL= ML + MMUON**2 - (gQ/2d0-g1Q/2d0)*(h2q**2-h1q**2)
      mslR= ME + MMUON**2 - g1Q/2d0*(h2q**2-h1q**2)
      Xl= AMUON-MUQ*tanbq
      Wl= DSQRT((mslL-mslR)**2 + 4d0*(Xl*MMUON)**2)
      MSMU1= 0.5d0*(mslL+mslR-Wl)
      MSMU2= 0.5d0*(mslL+mslR+Wl)

      nen= DSQRT((mslL-MSMU1)**2 + (Xl*MMUON)**2)
      IF(nen.EQ.0d0)THEN
       IF(mslL.LE.mslR)THEN
        CSMU= 1d0
       ELSE
        CSMU= 0d0
       ENDIF
      ELSE
       CSMU= -MMUON*Xl/nen
      ENDIF

      MSNM= ML + gQ/2d0*(h2q**2-h1q**2)

      IF(MSMU1.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MSMU1=1d2
!      WRITE(0,*)"MSMU1^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MSMU2.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MSMU2=1d2
!      WRITE(0,*)"MSMU2^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MSNM.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MSNM=1d2
!      WRITE(0,*)"MSNM^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MSMU1.GE.0d0)MSMU1=DSQRT(MSMU1)
      IF(MSMU1.GE.0d0)MSMU2=DSQRT(MSMU2)
      IF(MSNM.GE.0d0)MSNM=DSQRT(MSNM)

*   Selectron masses squared

      MLR=ME - g1Q/2d0*(h2q**2-h1q**2)
      MLL=ML + (-gQ/2d0+g1Q/2d0)*(h2q**2-h1q**2)
      MNL=ML + gQ/2d0*(h2q**2-h1q**2)

      IF(MLR.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MLR=1d2
!      WRITE(0,*)"MLR^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MLL.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MLL=1d2
!      WRITE(0,*)"MLL^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MNL.LE.1d2)THEN
       IFAIL=8
       IF(FLAG.EQ.0)MNL=1d2
!      WRITE(0,*)"MNL^2 < 0"
!      WRITE(0,*)""
!      WRITE(0,*)""
      ENDIF

      IF(MLR.GE.0d0)MLR=DSQRT(MLR)
      IF(MLL.GE.0d0)MLL=DSQRT(MLL)
      IF(MNL.GE.0d0)MNL=DSQRT(MNL)

      END

      DOUBLE PRECISION FUNCTION DMSQUARK(i,msq1,msq2,s2t,mq,mglu,q2)

*    with thanks to S. Kraml
*    msq1, msq2 are the squark masses squared

      IMPLICIT NONE
      DOUBLE PRECISION msq1,msq2,s2t,mq,mglu,msq,msqp
      DOUBLE PRECISION mq2,mglu2,sgn,dmg,dmsg,dmsq,q2
      DOUBLE PRECISION NMA0,NMB0,NMB1
      INTEGER i

      mglu=dabs(mglu)
      mq2 = mq*mq
      mglu2 = mglu*mglu
      IF (i.eq.1) THEN
        msq=msq1
        msqp=msq2
        sgn = -1d0
      ELSE
        msq=msq2
        msqp=msq1
        sgn = 1d0
      ENDIF

      dmg = -2d0*msq*( 2d0*NMB0(msq,0d0,msq,q2)-NMB1(msq,0d0,msq,q2) )
      dmsg = -4d0*( NMA0(mq2,q2) - msq*NMB1(msq,mglu2,mq2,q2) +
     .    (mglu2 + sgn*mglu*mq*s2t)*NMB0(msq,mglu2,mq2,q2) )
      dmsq = (1d0-s2t**2)*NMA0(msq,q2) + s2t**2*NMA0(msqp,q2)
      DMSQUARK = dmg + dmsg + dmsq

      END

      DOUBLE PRECISION FUNCTION INTEG(X,Y,Z)

* Function for DELMB:

      IMPLICIT NONE

      DOUBLE PRECISION X,Y,Z

       IF(DABS(X).EQ.DABS(Y) .AND. DABS(X).EQ.DABS(Z))THEN
        INTEG=.5d0/X
      ELSEIF(DABS(X).EQ.DABS(Y))THEN
        INTEG=(X**2-Z**2+Z**2*DLOG(Z**2/X**2))/(X**2-Z**2)**2
      ELSEIF(DABS(Y).EQ.DABS(Z))THEN
        INTEG=(Y**2-X**2+X**2*DLOG(X**2/Y**2))/(Y**2-X**2)**2
      ELSEIF(DABS(X).EQ.DABS(Z))THEN
        INTEG=(X**2-Y**2+Y**2*DLOG(Y**2/X**2))/(X**2-Y**2)**2
      ELSE
        INTEG=(X**2*Y**2*DLOG(X**2/Y**2)
     .   +Y**2*Z**2*DLOG(Y**2/Z**2)+Z**2*X**2*DLOG(Z**2/X**2))/
     .   ((X**2-Y**2)*(Y**2-Z**2)*(X**2-Z**2))
      ENDIF

      END

      SUBROUTINE DMSQUARK_YUK(IFAIL,DPI11ST,DPI22ST,DPI11SB,DPI22SB)

*******************************************************************
* Subroutine to compute the Yukawa contributions to the squark pole masses
* based on Pierce/Bagger et al., hep-ph/9606211 eqs. (D.47-49)
*
* NOTE: mu(Pierce/Bagger) = -MUQ(this code)!!!
*
* The conventions for the cosines/sines CSF/SSF (F=T, B) of the
* 3rd generation sfermion mixing angle are
* F1 = CSF*FL + SSF*FR, F2 = CSF*FR - SSF*FL
* where FR, FL are the left/right handed weak eigenstates, and
* F1, F2 are the physical states ordered in mass (M_F1 < M_F2)
*
* Comparing to (D.54) in hep-ph/9606211 one also finds
* F1 = CSF*FL + SSF*FR, F2 = CSF*FR - SSF*FL although F1, F2 are
* not ordered on mass, which is irrelevant
*
* Used trilin. couplings A_top, A_bottom:
* ATP, ABP at QSTSB from COMMON/RADCOR2/
* Used running stop/sbottom masses at QSTSB:
* RMST1/2, RMSB1/2 from COMMON/RADCOR/
*
*      SCOMP(1-3,1-3): Mixing angles: IF HB(I) are the weak eigenstates,
*        HB(I) = Re(H1), Re(H2), Re(S) where H1 == H_u, H2 == H_d
*        and HM(I) are the mass eigenstates,
*        the convention is HB(I) = SUM_(J=1,3) SCOMP(J,I)*HM(J)
*        which is equivalent to HM(I) = SUM_(J=1,3) SCOMP(I,J)*HB(J)
* -> SCOMP(K,1)) is the H_u-component of the mass eigenstate K
* -> SCOMP(K,2)) is the H_d-component of the mass eigenstate K
*
*      PCOMP(1-2,1-2): Mixing angles: IF AB(I) are the weak eigenstates,
*        AB(I) = Im(H1), Im(H2), Im(S) where A1 == A_u, A2 == A_d
*        and AM(I) are the mass eigenstates,
*        the convention is
*        AM(I) = PCOMP(I,1)*(COSBETA*AB(1)+SINBETA*AB(2))
*              + PCOMP(I,2)*AB(3)
* or, using
*       P(I,1)= PCOMP(I,1)*CB
*       P(I,2)= PCOMP(I,1)*SB
*       P(I,3)= PCOMP(I,2)
*       AM(I) = P(I,1)*AB(1)+P(I,2)*AB(2)) + P(I,3)*AB(3)

*******************************************************************

      IMPLICIT NONE

      INTEGER IFAIL,I,J,K,N

      DOUBLE PRECISION QSTSB,pi,coef
      DOUBLE PRECISION MUR,MUL,MDR,MDL,MLR,MLL,MNL
      DOUBLE PRECISION MST1,MST2,MSB1,MSB2,MSL1,MSL2,MSNT
      DOUBLE PRECISION CST,CSB,CSL,MSMU1,MSMU2,MSNM,CSMU
      DOUBLE PRECISION RMST1,RMST2,S2T,RMSB1,RMSB2,S2B,XT,XB
      DOUBLE PRECISION RMST(2),RMSB(2)
      DOUBLE PRECISION M1,M2,M3,MQ3,MU3,MD3,ML3,ME3,MQ,MU,MD,ML,ME
      DOUBLE PRECISION MQ3P,MU3P,MD3P,ATP,ABP
      DOUBLE PRECISION G1Q,G2Q,GQ,ALSQ
      DOUBLE PRECISION ZHU,ZHD,ZS,H1Q,H2Q,TANBQ
      DOUBLE PRECISION HTQ,HBQ,MTQ,MBQ
      DOUBLE PRECISION LQ,KQ,ALQ,AKQ,MUQ,NUQ
      DOUBLE PRECISION XI,HT2,HB2,MH1,MH2,L2,K2,HTAU2,AL,AK
      DOUBLE PRECISION SCOMP(3,3),PCOMP(2,2),CMASS,P(2,3)
      DOUBLE PRECISION H(3,3),A(2,2),MH(3),MA(2),VEC3(3,3),VEC2(2,2),EPS
      DOUBLE PRECISION XIF,XIS,MUP,MSP,M3H
      DOUBLE PRECISION RT,RB,FMT,FMB,AlshIFt,B,GMCOMB,SQR2

      DOUBLE PRECISION DPI11ST,DPI22ST,DPI11SB,DPI22SB
      DOUBLE PRECISION NMA0,NMB0,NMB1
      DOUBLE PRECISION A0T1,A0T2,A0B1,A0B2,A0H0(6),A0HC(2)
      DOUBLE PRECISION DNUST(6),DNDST(2),DNUSB(6),DNDSB(2)
      DOUBLE PRECISION B011HNTI(6,2),B022HNTI(6,2)
      DOUBLE PRECISION B011HCBI(2,2),B022HCBI(2,2)
      DOUBLE PRECISION B011HNBI(6,2),B022HNBI(6,2)
      DOUBLE PRECISION B011HCTI(2,2),B022HCTI(2,2)

      DOUBLE PRECISION SST,SSB,CB,SB
      DOUBLE PRECISION HRULUL,HRDLDL,HRURUR,HRDRDR,HRULUR
      DOUBLE PRECISION HRDLDR,HIULUR,HIDLDR,HGBULUR,HGBDLDR
      DOUBLE PRECISION LHNTLT(6,2),LHNTRT(6,2),LHNBLB(6,2),LHNBRB(6,2)
      DOUBLE PRECISION LHCTLBL(2),LHCTRBL(2),LHCTLBR(2),LHCTRBR(2)
      DOUBLE PRECISION LHCTLB(2,2),LHCTRB(2,2)
      DOUBLE PRECISION LHNBLT(6,2),LHNBRT(6,2),LHNTLB(6,2),LHNTRB(6,2)
      DOUBLE PRECISION LHCBLTL(2),LHCBRTL(2),LHCBLTR(2),LHCBRTR(2)
      DOUBLE PRECISION LHCBLT(2,2),LHCBRT(2,2)
      DOUBLE PRECISION MS,MC,MB,MBP,MT,MTAU,MMUON,MZ,MW

      COMMON/SFSPEC/MUR,MUL,MDR,MDL,MLR,MLL,MNL,
     .      MST1,MST2,MSB1,MSB2,MSL1,MSL2,MSNT,
     .      CST,CSB,CSL,MSMU1,MSMU2,MSNM,CSMU
      COMMON/RADCOR/RMST1,RMST2,S2T,RMSB1,RMSB2,S2B,XT,XB
      COMMON/RADCOR2/MQ3P,MU3P,MD3P,ATP,ABP
      COMMON/STSBSCALE/QSTSB
      COMMON/QGAUGE/G1Q,G2Q,GQ,ALSQ
      COMMON/QHIGGS/ZHU,ZHD,ZS,H1Q,H2Q,TANBQ
      COMMON/QQUARK/HTQ,HBQ,MTQ,MBQ
      COMMON/QPAR/LQ,KQ,ALQ,AKQ,MUQ,NUQ
      COMMON/QEXT/XIF,XIS,MUP,MSP,M3H
      COMMON/SMSPEC/MS,MC,MB,MBP,MT,MTAU,MMUON,MZ,MW

* DEFINITION OF COUPLINGS AS FUNCTIONS:

      HRULUL(K)= SQR2*(HTQ**2*H1Q*SCOMP(K,1))
      HRDLDL(K)= SQR2*(HBQ**2*H2Q*SCOMP(K,2))
      HRURUR(K)= SQR2*(HTQ**2*H1Q*SCOMP(K,1))
      HRDRDR(K)= SQR2*(HBQ**2*H2Q*SCOMP(K,2))
      HRULUR(K)= HTQ/SQR2*(ATP*SCOMP(K,1)-MUQ*SCOMP(K,2)
     . - LQ*H2Q*SCOMP(K,3))
      HRDLDR(K)= HBQ/SQR2*(-MUQ*SCOMP(K,1)+ABP*SCOMP(K,2)
     . - LQ*H1Q*SCOMP(K,3))
      HIULUR(K)= HTQ/SQR2*(ATP*P(K,1)+MUQ*P(K,2)
     . + LQ*H2Q*P(K,3))
      HIDLDR(K)= HBQ/SQR2*(MUQ*P(K,1)+ABP*P(K,2)
     . + LQ*H1Q*P(K,3))
* END: DEFINITION OF COUPLINGS AS FUNCTIONS

      EPS=1d-8
      PI=4d0*DATAN(1d0)
      COEF=1d0/(16d0*PI**2)
      SQR2= DSQRT(2d0)

      SST= DSQRT(1d0-CST**2)
      SSB= DSQRT(1d0-CSB**2)

*   Trig. functions of beta
      cb=1d0/DSQRT(1d0+tanbq**2)
      sb=tanbq*cb

********************************************
* First: Need an estimate of Higgs masses and mixings,
* tree level + top/stop-contr.; copied from mhiggs.f:
      rt= 3d0/2d0*COEF*htq**2
      rb= 3d0/2d0*COEF*hbq**2
       fmt= (rmst2*DLOG(rmst2/QSTSB)-rmst1*DLOG(rmst1/QSTSB))/
     .     (rmst2-rmst1)-1d0

       fmb= (rmsb2*DLOG(rmsb2/QSTSB)-rmsb1*DLOG(rmsb1/QSTSB))/
     .    (rmsb2-rmsb1)-1d0
      Alshift= ALQ+2d0*rt*Atp*fmt+2d0*rb*Abp*fmb
      B= Alshift+NUQ
      GMCOMB= LQ*XIF+MUP*MUQ+M3H

*   Tree level CP-even Higgs mass matrix
      H(1,1) = GQ*h1q**2 + (MUQ*B+GMCOMB)/tanbq
      H(1,1)= H(1,1) + rt*(4d0*mtq**2*DLOG(rmst1*rmst2/mtq**4))
      H(2,2) = GQ*h2q**2 + (MUQ*B+GMCOMB)*tanbq
      H(3,3) = LQ**2*(AlshIFt+MUP)*h1q*h2q/MUQ
     .   + NUQ*(AKQ+4d0*NUQ+3d0*MUP)-LQ/MUQ*(XIS+XIF*MUP)
      H(1,2) = (2d0*LQ**2-GQ)*h1q*h2q - MUQ*B-GMCOMB
      H(1,3) = LQ*(2d0*MUQ*h1q - (B+NUQ+MUP)*h2q)
      H(2,3) = LQ*(2d0*MUQ*h2q - (B+NUQ+MUP)*h1q)
*   Diagonalization
      CALL DIAGN(3,H,MH,VEC3,EPS)
      CALL SORTN(3,MH,VEC3)
      DO I= 1,3
       DO J= 1,3
        SCOMP(I,J)= VEC3(J,I)
       ENDDO
      ENDDO

*   CP-odd Higgs mass matrix including
      A(1,1)= (MUQ*B+GMCOMB)*(TANBQ+1D0/TANBQ)
      A(2,2)= (LQ**2*h1q*h2q/MUQ*(B+3d0*NUQ+MUP)-3d0*AKQ*NUQ
     .       -XIF*(4d0*KQ+LQ*MUP/MUQ)-2d0*MSP-MUP*NUQ-LQ*XIS/MUQ)
      A(1,2)= LQ*(ALQ+2d0*rt*Atp*fmt+2d0*rb*Abp*fmb)
     .       *DSQRT(h1q**2+h2q**2)
*   Diagonalization
      CALL DIAGN(2,A,MA,VEC2,EPS)
      CALL SORTN(2,MA,VEC2)
      DO I= 1,2
       DO J= 1,2
        PCOMP(I,J)= VEC2(J,I)
       ENDDO
      ENDDO
*   CP odd mixing angles
      DO I=1,2
       P(I,1)= PCOMP(I,1)*CB
       P(I,2)= PCOMP(I,1)*SB
       P(I,3)= PCOMP(I,2)
      ENDDO

*   Charged Higgs mass
      CMASS=((G2Q/2d0-LQ**2)*h1q*h2q+MUQ*B+GMCOMB)*(TANBQ+1D0/TANBQ)

      DO I=1,3
       MH(I)=MAX(MH(I),10d0)
      ENDDO
      DO I=1,2
       MA(I)=MAX(MA(I),10d0)
      ENDDO
      CMASS=MAX(CMASS,10d0)

* END estimate of Higgs masses
********************************************

      RMST(1)=MAX(RMST1,10d0)
      RMST(2)=MAX(RMST2,10d0)
      RMSB(1)=MAX(RMSB1,10d0)
      RMSB(2)=MAX(RMSB2,10d0)

      A0T1=NMA0(RMST(1),QSTSB)
      A0T2=NMA0(RMST(2),QSTSB)
      A0B1=NMA0(RMSB(1),QSTSB)
      A0B2=NMA0(RMSB(2),QSTSB)

* DEFINITIONS OF COUPLINGS AND PASSARINO-VELTMAN-FUNCTIONS
* defined in subfun.f
* CP-even neutral Higgses: J: MH(1-3)
      DO J=1,3
          A0H0(J)=NMA0(MH(J),QSTSB)
          DNUST(J)=SCOMP(J,1)**2
          DNUSB(J)=SCOMP(J,2)**2
* Couplings lambda_Higgs(J)_STOP(L,R)_STOP(1,2):
          LHNTLT(J,1)=CST*HRULUL(J)+SST*HRULUR(J)
          LHNTLT(J,2)=-SST*HRULUL(J)+CST*HRULUR(J)
          LHNTRT(J,1)=CST*HRULUR(J)+SST*HRURUR(J)
          LHNTRT(J,2)=-SST*HRULUR(J)+CST*HRURUR(J)
* Couplings lambda_Higgs(J)_SBOT(L,R)_SBOT(1,2):
          LHNBLB(J,1)=CSB*HRDLDL(J)+SSB*HRDLDR(J)
          LHNBLB(J,2)=-SSB*HRDLDL(J)+CSB*HRDLDR(J)
          LHNBRB(J,1)=CSB*HRDLDR(J)+SSB*HRDRDR(J)
          LHNBRB(J,2)=-SSB*HRDLDR(J)+CSB*HRDRDR(J)
* I: stop1/stop2 or sbot1/sbot2
        DO I=1,2
          B011HNTI(J,I)=NMB0(RMST(1),MH(J),RMST(I),QSTSB)
          B022HNTI(J,I)=NMB0(RMST(2),MH(J),RMST(I),QSTSB)
          B011HNBI(J,I)=NMB0(RMSB(1),MH(J),RMSB(I),QSTSB)
          B022HNBI(J,I)=NMB0(RMSB(2),MH(J),RMSB(I),QSTSB)
        ENDDO
      ENDDO
* END CP-even neutral Higgses
* CP-odd neutral Higgses:  J: MA(1-2)
      DO J=1,2
          A0H0(3+J)=NMA0(MA(J),QSTSB)
          DNUST(3+J)=P(J,1)**2
          DNUSB(3+J)=P(J,2)**2
* Couplings lambda_Higgs(J)_STOP(L,R)_STOP(1,2):
          LHNTLT(3+J,1)=SST*HIULUR(J)
          LHNTLT(3+J,2)=CST*HIULUR(J)
          LHNTRT(3+J,1)=CST*HIULUR(J)
          LHNTRT(3+J,2)=-SST*HIULUR(J)
* Couplings lambda_Higgs(J)_SBOT(L,R)_SBOT(1,2):
          LHNBLB(3+J,1)=SSB*HIDLDR(J)
          LHNBLB(3+J,2)=CSB*HIDLDR(J)
          LHNBRB(3+J,1)=CSB*HIDLDR(J)
          LHNBRB(3+J,2)=-SSB*HIDLDR(J)
* I: stop1/stop2 or sbot1/sbot2
        DO I=1,2
          B011HNTI(3+J,I)=NMB0(RMST(1),MA(J),RMST(I),QSTSB)
          B022HNTI(3+J,I)=NMB0(RMST(2),MA(J),RMST(I),QSTSB)
          B011HNBI(3+J,I)=NMB0(RMSB(1),MA(J),RMSB(I),QSTSB)
          B022HNBI(3+J,I)=NMB0(RMSB(2),MA(J),RMSB(I),QSTSB)
        ENDDO
      ENDDO
* END CP-odd neutral Higgses
* Couplings of the neutral Goldstone:
      HGBULUR=HTQ/SQR2*(ATP*SB-MUQ*CB)
      HGBDLDR= HBQ/SQR2*(MUQ*SB-ABP*CB)
* J=3: Goldstone Boson
          A0H0(3+3)=NMA0(MZ**2,QSTSB)
          DNUST(3+3)=SB**2
          DNUSB(3+3)=CB**2
* Couplings lambda_Higgs(J)_STOP(L,R)_STOP(1,2):
          LHNTLT(3+3,1)=SST*HGBULUR
          LHNTLT(3+3,2)=CST*HGBULUR
          LHNTRT(3+3,1)=CST*HGBULUR
          LHNTRT(3+3,2)=-SST*HGBULUR
* Couplings lambda_Higgs(J)_SBOT(L,R)_SBOT(1,2):
          LHNBLB(3+3,1)=SSB*HGBDLDR
          LHNBLB(3+3,2)=CSB*HGBDLDR
          LHNBRB(3+3,1)=CSB*HGBDLDR
          LHNBRB(3+3,2)=-SSB*HGBDLDR
* I: stop1/stop2 or sbot1/sbot2
        DO I=1,2
          B011HNTI(3+3,I)=NMB0(RMST(1),MZ**2,RMST(I),QSTSB)
          B022HNTI(3+3,I)=NMB0(RMST(2),MZ**2,RMST(I),QSTSB)
          B011HNBI(3+3,I)=NMB0(RMSB(1),MZ**2,RMSB(I),QSTSB)
          B022HNBI(3+3,I)=NMB0(RMSB(2),MZ**2,RMSB(I),QSTSB)
        ENDDO
* END Goldstone Boson
* Charged Higgs:
        A0HC(1) = NMA0(CMASS,QSTSB)
        DNDST(1)=CB**2
        DNDSB(1)=SB**2
* Couplings lambda_CHiggs_STOP(L,R)_SBOT(L/R):
        LHCTLBL(1)=-HTQ**2*H1Q*CB-HBQ**2*H2Q*SB
        LHCTRBR(1)=-HTQ*HBQ*(H2Q*CB+H1Q*SB)
        LHCTLBR(1)=-HBQ*(MUQ*CB+ABP*SB)
        LHCTRBL(1)=-HTQ*(MUQ*SB+ATP*CB)
* Couplings lambda_CHiggs_SBOT(L,R)_STOP(L/R):
        LHCBLTL(1)=-HTQ**2*H1Q*CB-HBQ**2*H2Q*SB
        LHCBRTR(1)=-HTQ*HBQ*(H2Q*CB+H1Q*SB)
        LHCBLTR(1)=-HTQ*(MUQ*SB+ATP*CB)
        LHCBRTL(1)=-HBQ*(MUQ*CB+ABP*SB)
* LHCTLB(i,j): i=1: charged Higgs; j: sbottom1/2
        LHCTLB(1,1)= CSB*LHCTLBL(1)+SSB*LHCTLBR(1)
        LHCTLB(1,2)=-SSB*LHCTLBL(1)+CSB*LHCTLBR(1)
        LHCTRB(1,1)= CSB*LHCTRBL(1)+SSB*LHCTRBR(1)
        LHCTRB(1,2)=-SSB*LHCTRBL(1)+CSB*LHCTRBR(1)
* LHCBLT(i,j): i=1: charged Higgs; j: stop1/2
        LHCBLT(1,1)= CST*LHCBLTL(1)+SST*LHCBLTR(1)
        LHCBLT(1,2)=-SST*LHCBLTL(1)+CST*LHCBLTR(1)
        LHCBRT(1,1)= CST*LHCBRTL(1)+SST*LHCBRTR(1)
        LHCBRT(1,2)=-SST*LHCBRTL(1)+CST*LHCBRTR(1)
* I: sbot1/sbot2 or stop1/stop2
      DO I=1,2
        B011HCBI(1,I)=NMB0(RMST(1),CMASS,RMSB(I),QSTSB)
        B022HCBI(1,I)=NMB0(RMST(2),CMASS,RMSB(I),QSTSB)
        B011HCTI(1,I)=NMB0(RMSB(1),CMASS,RMST(I),QSTSB)
        B022HCTI(1,I)=NMB0(RMSB(2),CMASS,RMST(I),QSTSB)
      ENDDO
* Charged Goldstone:
        A0HC(2) = NMA0(MW**2,QSTSB)
        DNDST(2)=SB**2
        DNDSB(2)=CB**2
        LHCTLBL(2)=-HTQ**2*H1Q*SB+HBQ**2*H2Q*CB
        LHCTRBR(2)=0D0
        LHCTLBR(2)=HBQ*(-MUQ*SB+ABP*CB)
        LHCTRBL(2)=HTQ*(MUQ*CB-ATP*SB)
        LHCBLTL(2)=-HTQ**2*H1Q*SB+HBQ**2*H2Q*CB
        LHCBRTR(2)=0D0
        LHCBLTR(2)=HTQ*( MUQ*CB-ATP*SB)
        LHCBRTL(2)=HBQ*(-MUQ*SB+ABP*CB)
* LHCTLB(i,j): i=2: charged Goldstone; j: sbottom1/2
        LHCTLB(2,1)=CSB*LHCTLBL(2)+SSB*LHCTLBR(2)
        LHCTLB(2,2)=-SSB*LHCTLBL(2)+CSB*LHCTLBR(2)
        LHCTRB(2,1)=CSB*LHCTRBL(2)+SSB*LHCTRBR(2)
        LHCTRB(2,2)=-SSB*LHCTRBL(2)+CSB*LHCTRBR(2)
* LHCBLT(i,j): i=2: charged Goldstone; j: stop1/2
        LHCBLT(2,1)= CST*LHCBLTL(2)+SST*LHCBLTR(2)
        LHCBLT(2,2)=-SST*LHCBLTL(2)+CST*LHCBLTR(2)
        LHCBRT(2,1)= CST*LHCBRTL(2)+SST*LHCBRTR(2)
        LHCBRT(2,2)=-SST*LHCBRTL(2)+CST*LHCBRTR(2)
* I: sbot1/sbot2 or stop1/stop2
      DO I=1,2
        B011HCBI(2,I)=NMB0(RMST(1),MW**2,RMSB(I),QSTSB)
        B022HCBI(2,I)=NMB0(RMST(2),MW**2,RMSB(I),QSTSB)
        B011HCTI(2,I)=NMB0(RMSB(1),MW**2,RMST(I),QSTSB)
        B022HCTI(2,I)=NMB0(RMSB(2),MW**2,RMST(I),QSTSB)
      ENDDO
* END Charged Goldstone

* CONTRIBUTIONS TO DPI_11ST:
      DPI11ST=HTQ**2*(2D0*(CST*SST)**2*A0T1+(CST**4+SST**4)*A0T2)
     .       +((HBQ**2*(CST*SSB)**2+HTQ**2*(SST*CSB)**2)*A0B1
     .        +(HBQ**2*(CST*CSB)**2+HTQ**2*(SST*SSB)**2)*A0B2)
     .       +6d0*HTQ**2*CST**2*SST**2*(A0T1-A0T2)
*  N=1,2,3: CP-EVEN; I=4,5: CP-ODD; I=6: GOLDSTONE
      DO N=1,6
        DPI11ST=DPI11ST+HTQ**2/2D0*DNUST(N)*A0H0(N)

*   I: STOP1/STOP2
        DO I=1,2
        If(N.le.3)then
          DPI11ST=DPI11ST+(CST*LHNTLT(N,I)+SST*LHNTRT(N,I))**2
     .                   *B011HNTI(N,I)
        else
          DPI11ST=DPI11ST+(CST*LHNTLT(N,I)-SST*LHNTRT(N,I))**2
     .                   *B011HNTI(N,I)
        endif
        ENDDO
      ENDDO
*  N=1: CHARGED HIGGS; N=2: CHARGED GOLDSTONE
      DO N=1,2
        DPI11ST=DPI11ST
     . +(HBQ**2*CST**2*DNDST(N)+HTQ**2*SST**2*DNUST(N))*A0HC(N)
*   I: SBOT1/SBOT2
        DO I=1,2
          DPI11ST=DPI11ST+(CST*LHCTLB(N,I)+SST*LHCTRB(N,I))**2
     .                   *B011HCBI(N,I)
        ENDDO
      ENDDO

        DPI11ST=DPI11ST+HTQ**2
     . *((RMST(1)-muq**2-mtq**2)*NMB0(RMST(1),muq**2,mtq**2,QSTSB)
     .                       -NMA0(muq**2,QSTSB)-NMA0(mtq**2,QSTSB))
     . +(HBQ**2*CST**2+HTQ**2*SST**2)
     . *((RMST(1)-muq**2-mbq**2)*NMB0(RMST(1),muq**2,mbq**2,QSTSB)
     .                       -NMA0(muq**2,QSTSB)-NMA0(mbq**2,QSTSB))
     . -4d0*CST*SST
     .     *HTQ*HBQ*muq*mbq*NMB0(RMST(1),muq**2,mbq**2,QSTSB)
      DPI11ST=DPI11ST*COEF

* CONTRIBUTIONS TO DPI_11SB:
      DPI11SB=HBQ**2*(2D0*(CSB*SSB)**2*A0B1+(CSB**4+SSB**4)*A0B2)
     .       +((HTQ**2*(CSB*SST)**2+HBQ**2*(SSB*CST)**2)*A0T1
     .        +(HTQ**2*(CSB*CST)**2+HBQ**2*(SSB*SST)**2)*A0T2)
     .       +6d0*HBQ**2*CSB**2*SSB**2*(A0B1-A0B2)
*  N=1,2,3: CP-EVEN; I=4,5: CP-ODD; I=6: GOLDSTONE
      DO N=1,6
        DPI11SB=DPI11SB+HBQ**2/2D0*DNUSB(N)*A0H0(N)
*   I: SBOT1/SBOT2
        DO I=1,2
        If(N.le.3)then
          DPI11SB=DPI11SB+(CSB*LHNBLB(N,I)+SSB*LHNBRB(N,I))**2
     .                   *B011HNBI(N,I)
        else
          DPI11SB=DPI11SB+(CSB*LHNBLB(N,I)-SSB*LHNBRB(N,I))**2
     .                   *B011HNBI(N,I)
        endif
        ENDDO
      ENDDO
*  N=1: CHARGED HIGGS; N=2: CHARGED GOLDSTONE
      DO N=1,2
        DPI11SB=DPI11SB
     .      +(HTQ**2*CSB**2*DNUSB(N)+HBQ**2*SSB**2*DNDSB(N))*A0HC(N)
*   I: STOP1/STOP2
        DO I=1,2
          DPI11SB=DPI11SB+(CSB*LHCBLT(N,I)+SSB*LHCBRT(N,I))**2
     .                   *B011HCTI(N,I)
        ENDDO
      ENDDO

      DPI11SB=DPI11SB+HBQ**2
     . *((RMSB(1)-muq**2-mbq**2)*NMB0(RMSB(1),muq**2,mbq**2,QSTSB)
     .                       -NMA0(muq**2,QSTSB)-NMA0(mbq**2,QSTSB))
     . +(HTQ**2*CSB**2+HBQ**2*SSB**2)
     . *((RMSB(1)-muq**2-mtq**2)*NMB0(RMSB(1),muq**2,mtq**2,QSTSB)
     .                       -NMA0(muq**2,QSTSB)-NMA0(mtq**2,QSTSB))
     . -4d0*CSB*SSB
     .     *HTQ*HBQ*muq*mtq*NMB0(RMSB(1),muq**2,mtq**2,QSTSB)
      DPI11SB=DPI11SB*COEF

* CONTRIBUTIONS TO DPI_22ST:
      DPI22ST=HTQ**2*(2D0*(CST*SST)**2*A0T2+(CST**4+SST**4)*A0T1)
     .       +((HTQ**2*(CST*SSB)**2+HBQ**2*(SST*CSB)**2)*A0B2
     .        +(HTQ**2*(CST*CSB)**2+HBQ**2*(SST*SSB)**2)*A0B1)
     .       +6d0*HTQ**2*CST**2*SST**2*(A0T2-A0T1)
*  N=1,2,3: CP-EVEN; I=4,5: CP-ODD; I=6: GOLDSTONE
      DO N=1,6
        DPI22ST=DPI22ST+HTQ**2/2D0*DNUST(N)*A0H0(N)
*   I: STOP1/STOP2
        DO I=1,2
        if(N.le.3)then
          DPI22ST=DPI22ST+(-SST*LHNTLT(N,I)+CST*LHNTRT(N,I))**2
     .                   *B022HNTI(N,I)
        else
          DPI22ST=DPI22ST+(-SST*LHNTLT(N,I)-CST*LHNTRT(N,I))**2
     .                   *B022HNTI(N,I)
        endif
        ENDDO
      ENDDO
*  N=1: CHARGED HIGGS; N=2: CHARGED GOLDSTONE
      DO N=1,2
        DPI22ST=DPI22ST
     .   +(HBQ**2*SST**2*DNDST(N)+HTQ**2*CST**2*DNUST(N))*A0HC(N)
*   I: SBOT1/SBOT2
        DO I=1,2
          DPI22ST=DPI22ST+(-SST*LHCTLB(N,I)+CST*LHCTRB(N,I))**2
     .                   *B022HCBI(N,I)
        ENDDO
      ENDDO

        DPI22ST=DPI22ST+HTQ**2
     . *((RMST(2)-muq**2-mtq**2)*NMB0(RMST(2),muq**2,mtq**2,QSTSB)
     .                       -NMA0(muq**2,QSTSB)-NMA0(mtq**2,QSTSB))
     . +(HBQ**2*SST**2+HTQ**2*CST**2)
     . *((RMST(2)-muq**2-mbq**2)*NMB0(RMST(2),muq**2,mbq**2,QSTSB)
     .                       -NMA0(muq**2,QSTSB)-NMA0(mbq**2,QSTSB))
     . +4d0*CST*SST
     .     *HTQ*HBQ*muq*mbq*NMB0(RMST(2),muq**2,mbq**2,QSTSB)
      DPI22ST=DPI22ST*COEF

* CONTRIBUTIONS TO DPI_22SB:
      DPI22SB=HBQ**2*(2D0*(CSB*SSB)**2*A0B2+(CSB**4+SSB**4)*A0B1)
     .       +((HBQ**2*(CSB*SST)**2+HTQ**2*(SSB*CST)**2)*A0T2
     .        +(HBQ**2*(CSB*CST)**2+HTQ**2*(SSB*SST)**2)*A0T1)
     .       +6d0*HBQ**2*CSB**2*SSB**2*(A0B2-A0B1)
*  N=1,2,3: CP-EVEN; I=4,5: CP-ODD; I=6: GOLDSTONE
      DO N=1,6
        DPI22SB=DPI22SB+HBQ**2/2D0*DNUSB(N)*A0H0(N)
*   I: SBOT1/SBOT2
        DO I=1,2
        if(N.le.3)then
          DPI22SB=DPI22SB+(-SSB*LHNBLB(N,I)+CSB*LHNBRB(N,I))**2
     .                   *B022HNBI(N,I)
        else
          DPI22SB=DPI22SB+(-SSB*LHNBLB(N,I)-CSB*LHNBRB(N,I))**2
     .                   *B022HNBI(N,I)
        endif
        ENDDO
      ENDDO
*  N=1: CHARGED HIGGS; N=2: CHARGED GOLDSTONE
      DO N=1,2
        DPI22SB=DPI22SB
     .    +(HTQ**2*SSB**2*DNUSB(N)+HBQ**2*CSB**2*DNDSB(N))*A0HC(N)
*   I: STOP1/STOP2
        DO I=1,2
          DPI22SB=DPI22SB+(-SSB*LHCBLT(N,I)+CSB*LHCBRT(N,I))**2
     .                   *B022HCTI(N,I)
        ENDDO
      ENDDO

      DPI22SB=DPI22SB+HBQ**2
     . *((RMSB(2)-muq**2-mbq**2)*NMB0(RMSB(2),muq**2,mbq**2,QSTSB)
     .                       -NMA0(muq**2,QSTSB)-NMA0(mbq**2,QSTSB))
     . +(HTQ**2*SSB**2+HBQ**2*CSB**2)
     . *((RMSB(2)-muq**2-mtq**2)*NMB0(RMSB(2),muq**2,mtq**2,QSTSB)
     .                       -NMA0(muq**2,QSTSB)-NMA0(mtq**2,QSTSB))
     . +4d0*CSB*SSB
     .     *HTQ*HBQ*muq*mtq*NMB0(RMSB(2),muq**2,mtq**2,QSTSB)
      DPI22SB=DPI22SB*COEF

      END