text:
      SUBROUTINE CHECKMIN_CPV(PAR,PROB)

**********************************************************************
c The effective potential (up to a constant term) is evaluated at the
c electroweak minimum, VEW. It is then compared to the value of the
c potential at several other points:
c  * VZ3: for <S>=<Hu>=<Hd>=0
c  * Vgauge: for <Hu>=<Hd>=0, <S>=mu/l
c  * Vvu: for <S>=<Hd>=0, <Hu>=vu
c  * Vvd: for <S>=<Hu>=0, <Hd>=vd
c  * Vphi: for <S>=mu/l.exp(i phiS), <Hu.Hd>=vu.vd.exp(i phiH)
c where phiS and phiH are scanned over on a Nphi*Nphi grid.
c Should one of these minima prove deeper than VEW, PROB(28) is.
c switched on
c The quadratic terms of the potential mHu2 and mHd2 are computed and
c it is checked whether they lie too far from the SUSY scale, in which
c case PROB(29) neq 0.
**********************************************************************

      IMPLICIT NONE

      INTEGER I,J,Nphi,CFLAG(5)

      DOUBLE PRECISION Pi,aux,Taux,Daux,Ytau
      DOUBLE PRECISION MW2,MZ2,Fsf1
      DOUBLE PRECISION VEW,VZ3,Vgauge,Vvu,Vvd,Vphi
      DOUBLE PRECISION phiS,phiH
      DOUBLE PRECISION MSTaux1,MSTaux2,MSBaux1,MSBaux2,MSLaux1,MSLaux2,
     . MSNTaux,MSUaux1,MSUaux2,MSDaux1,MSDaux2,MSEaux1,MSEaux2,MSNEaux
      DOUBLE PRECISION MHD1aux,MHD2aux,MS1aux,MS2aux
      DOUBLE PRECISION MH2(6,6),VALPH(6),VECPH(6,6),Mch1aux,Mch2aux
      DOUBLE PRECISION MN2(10,10),VALPN(10),VECPN(10,10)
      DOUBLE PRECISION Aroot,Broot,Croot,Droot,Proot,Qroot,Delt,root1,
     . root2,root3

      DOUBLE PRECISION QSTSB
      DOUBLE PRECISION Q2
      DOUBLE PRECISION G1Q,G2Q,GQ,ALSQ
      DOUBLE PRECISION tanb,cosb,sinb,vu,vd
      DOUBLE PRECISION l,k,Alcos1,Akcos2,muq,nuq
      DOUBLE PRECISION ZHU,ZHD,ZS,vuq,vdq,TANBQ
      DOUBLE PRECISION XIFQ,XISQ,MUPQ,MSPQ,M3HQ
      DOUBLE PRECISION mt,mb,mtau,mmu,mel,MS,MC,MBP,MPI,MSTRANGE
      DOUBLE PRECISION Ytq,Ybq,MTOPQ,MBOTQ
      DOUBLE PRECISION MST2(2),UT(2,2,2),MSB2(2),UB(2,2,2),MSL2(2),
     . UTAU(2,2,2),MSNT2
      DOUBLE PRECISION MSU2(2),MSD2(2),MSE2(2),MSNE2,MSMU2(2),
     . UMU(2,2,2)
      DOUBLE PRECISION MSQ3,MSU3,MSD3,ATP,ABP
      DOUBLE PRECISION MSQ1,MSU1,MSD1
      DOUBLE PRECISION MSL3,MSE3,MSL1,MSE1,ATAU,AMU
      DOUBLE PRECISION MCH2(2),U(2,2,2),V(2,2,2)
      DOUBLE PRECISION MNEU(5),NEU(5,5,2)
      DOUBLE PRECISION mur,M1r,M2r,msi
      DOUBLE PRECISION MHC,XC(2,2),MH0(5),XH(5,5),MA2
      DOUBLE PRECISION PAR(*),PROB(*)
      DOUBLE PRECISION MHuS,MHdS,MSS
      DOUBLE PRECISION IAL,IAK,IXIS
      DOUBLE PRECISION phiF,phiP,phi3,phiSs,phiSP,phi3q,phiSq,phiSPq,
     .              mupsi,ks2si
      DOUBLE PRECISION phi01,phi02,phi0,phiM1,phiM2,phiM3,
     .              phiAT,phiAB,phiATAU,phiAC,phiAS,phiAMU
      DOUBLE PRECISION XIF,XIS,MUP,MSP,M3H
      DOUBLE PRECISION RXIS,RXIF,IAlQ2,IAkQ2,IXISQ2,MHuQ2,MHdQ2,MSQ2
      DOUBLE PRECISION DDCOS,DDSIN

      COMMON/PHASES/phi01,phi02,phi0,phiM1,phiM2,phiM3,
     .              phiAT,phiAB,phiATAU,phiAC,phiAS,phiAMU
      COMMON/STSBSCALE/QSTSB
      COMMON/RENSCALE/Q2
      COMMON/QGAUGE/G1Q,G2Q,GQ,ALSQ
      COMMON/TBPAR/tanb,cosb,sinb,vu,vd
      COMMON/QPAR/l,k,Alcos1,Akcos2,muq,NUQ
      COMMON/QHIGGS/ZHU,ZHD,ZS,vuq,vdq,TANBQ
      COMMON/QEXT/XIFQ,XISQ,MUPQ,MSPQ,M3HQ
      COMMON/SMFERM/mt,mb,mtau,mmu,mel,MS,MC,MBP,MPI,MSTRANGE
      COMMON/QQUARK/Ytq,Ybq,MTOPQ,MBOTQ
      COMMON/SFERM3SPEC/MST2,UT,MSB2,UB,MSL2,UTAU,MSNT2
      COMMON/SFERM1SPEC/MSU2,MSD2,MSE2,MSNE2,MSMU2,UMU
      COMMON/RADCOR2/MSQ3,MSU3,MSD3,ATP,ABP
      COMMON/SQUPAR/MSQ1,MSU1,MSD1
      COMMON/SLEPPAR/MSL3,MSE3,MSL1,MSE1,ATAU,AMU
      COMMON/CHASPEC/MCH2,U,V
      COMMON/NEUSPEC/MNEU,NEU
      COMMON/GAUGINOPAR/mur,M1r,M2r,msi
      COMMON/HISPEC/MHC,XC,MH0,XH,MA2
      COMMON/MH2TREE/MHuS,MHdS,MSS
      COMMON/IMALAK/IAL,IAK,IXIS
      COMMON/Z3VAUX/phiF,phiP,phi3,phiSs,phiSP,phi3q,phiSq,phiSPq,
     .              mupsi,ks2si
      COMMON/SUSYEXT/XIF,XIS,MUP,MSP,M3H
      COMMON/REXI/RXIS,RXIF,IAlQ2,IAkQ2,IXISQ2,MHuQ2,MHdQ2,MSQ2
      COMMON/CFLAG/CFLAG

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

      Pi=4d0*datan(1d0)

c      Running to Q2

      Ytau=mtau/vdq

      IAlQ2=IAl-(4d0*l**2*IAl
     .    +2d0*k**2*(IAk*DDCOS(phi0)+Akcos2*DDSIN(phi0))
     .    +3d0*Ytq**2*PAR(12)*DDSIN(phiAt+phi01)
     .    +3d0*Ybq**2*PAR(13)*DDSIN(phiAb+phi01)
     .    +Ytau**2*PAR(14)*DDSIN(phiAtau+phi01)
     .    +g1q*PAR(20)*DDSIN(phiM1+phi01)
     .    +3d0*g2q*PAR(21)*DDSIN(phiM2+phi01))*DLOG(QSTSB/Q2)/16d0/Pi**2

	If(k.ne.0d0)then
      IAkQ2=IAk-6d0*(l**2*(IAl*DDCOS(phi0)-PAR(5)*DDSIN(phi0))
     .     +k**2*IAk)*DLOG(QSTSB/Q2)/16d0/Pi**2

      IXISQ2=XIS*DDSIN(PhiSs)
	else
      IAkQ2=0d0

      IXISQ2=IXIS
     .        *(1d0-(L**2+K**2)/16d0/Pi**2*DLOG(QSTSB/Q2))
     .  -(2d0*XIF*(L**2*(IAl*DDCOS(phiF-phi01)+PAR(5)*DDSIN(phiF-phi01))
     .           +K**2*(IAk*DDCOS(phiF-phi02)+PAR(6)*DDSIN(phiF-phi02)))
     .   +2d0*L*M3H*(IAl*DDCOS(phi3)-PAR(5)*DDSIN(phi3)
     .                             +MUP*DDSIN(phi3+phiP-phi01))
     .   +K*MSP*(IAk*DDCOS(phiSP)-PAR(6)*DDSIN(phiSP)
     .                             +MUP*DDSIN(phip+phiSP-phi02))
     .   +2d0*k*mup*MSS*DDSIN(phi02-phiP))/16d0/Pi**2*DLOG(QSTSB/Q2)
	endif

      MHuQ2=MHuS-(3d0*Ytq**2*(MSQ3+MSU3+MHuS+ATP**2)
     .           +L**2*(MHuS+MHdS+MSS+Alcos1**2+IAl**2)
     .           -g1q*PAR(20)**2-3*g2q*PAR(21)**2
     .           +g1q/2d0*(2d0*MSQ1+MSQ3-4d0*MSU1-2d0*MSU3+2d0*MSD1+MSD3
     .                    -2d0*MSL1-MSL3+2d0*MSE1+MSE3+MHuS-MHdS))
     .           /16d0/Pi**2*DLOG(QSTSB/Q2)

      MHdQ2=MHdS-(3d0*Ybq**2*(MSQ3+MSD3+MHuS+ABP**2)
     .           +Ytau**2*(MSL3+MSE3+MHdS+ATAU**2)
     .           +L**2*(MHuS+MHdS+MSS+Alcos1**2+IAl**2)
     .           -g1q*PAR(20)**2-3*g2q*PAR(21)**2
     .           +g1q/2d0*(2d0*MSQ1+MSQ3-4d0*MSU1-2d0*MSU3+2d0*MSD1+MSD3
     .                    -2d0*MSL1-MSL3+2d0*MSE1+MSE3+MHuS-MHdS))
     .           /16d0/Pi**2*DLOG(QSTSB/Q2)

      MSQ2=MSS-(2d0*L**2*(MHuS+MHdS+MSS+Alcos1**2+IAl**2)
     .         +2d0*K**2*(3d0*MSS+Akcos2**2+IAk**2))
     .           /16d0/Pi**2*DLOG(QSTSB/Q2)

      IF(CFLAG(1).EQ.0)RETURN

c      Preliminaries

      Nphi=10
      MW2=g2q*(vuq**2+vdq**2)/2d0
      MZ2=(g1q+g2q)*(vuq**2+vdq**2)/2d0

c      Electroweak minimum

      VEW=MHuS*vuq**2+MHdS*vdq**2+MSS*(muq/l)**2-2d0*muq*vuq*vdq*Alcos1
     .    +2d0/3d0*k*Akcos2*(muq/l)**3+muq**2*(vuq**2+vdq**2)
     .    +(l*vuq*vdq)**2-2d0*k/l*muq**2*DDCOS(phi0)*vuq*vdq
     .    +k**2*(muq/l)**4+(G1Q+G2Q)/8d0*(vuq**2-vdq**2)**2
     .    -2d0*(M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-phiF)
     .                        +MUPQ*muq*DDCOS(phi01-phiP))*vuq*vdq
     .    +2d0*(XISQ*DDCOS(phiSq)+XIFQ*MUPQ*DDCOS(phiP-phiF))*muq/l
     .    +2d0*(MSPQ/2d0*DDCOS(phiSPq)+k*XIFQ*DDCOS(phi02-phiF))
     .     *(muq/l)**2+2d0*k*MUPQ*DDCOS(phi02-phiP)*(muq/l)**3
     .    +(MUPQ*muq/l)**2
     .    -(3d0*mtopq**4*(dlog(mtopq**2/QSTSB)-1.5d0)
     .     +3d0*mbotq**4*(dlog(mbotq**2/QSTSB)-1.5d0)
     .     +(mtau*vdq/vd)**4*(dlog(mtau**2/QSTSB)-1.5d0))/16d0/Pi**2
     .    +3d0*(2d0*MW2**2*(dlog(MW2/QSTSB)-1.5d0)
     .                +MZ2*(dlog(MZ2/QSTSB)-1.5d0))/64d0/Pi**2
     .    +(3d0*MST2(1)**2*(dlog(MST2(1)/QSTSB)-1.5d0)
     .     +3d0*MST2(2)**2*(dlog(MST2(2)/QSTSB)-1.5d0)
     .     +3d0*MSB2(1)**2*(dlog(MSB2(1)/QSTSB)-1.5d0)
     .     +3d0*MSB2(2)**2*(dlog(MSB2(2)/QSTSB)-1.5d0)
     .     +MSL2(1)**2*(dlog(MSL2(1)/QSTSB)-1.5d0)
     .     +MSL2(2)**2*(dlog(MSL2(2)/QSTSB)-1.5d0)
     .     +MSNT2**2*(dlog(MSNT2/QSTSB)-1.5d0)
     .     +6d0*MSU2(1)**2*(dlog(MSU2(1)/QSTSB)-1.5d0)
     .     +6d0*MSU2(2)**2*(dlog(MSU2(2)/QSTSB)-1.5d0)
     .     +6d0*MSD2(1)**2*(dlog(MSD2(1)/QSTSB)-1.5d0)
     .     +6d0*MSD2(2)**2*(dlog(MSD2(2)/QSTSB)-1.5d0)
     .     +2d0*MSE2(1)**2*(dlog(MSE2(1)/QSTSB)-1.5d0)
     .     +2d0*MSE2(2)**2*(dlog(MSE2(2)/QSTSB)-1.5d0)
     .     +2d0*MSNE2**2*(dlog(MSNE2/QSTSB)-1.5d0))/32d0/Pi**2
     .    +3d0/512d0/Pi**4*Ytq**4*((dlog(QSTSB/mtopq**2))**2
     .      *(64d0*Pi*ALSQ+4d0/3d0*g1q-3d0*sinb**2*Ytq**2
     .       +3d0*cosb**2*Ybq**2)+
     .       ((dlog(MA2/mtopq**2))**2-(dlog(QSTSB/mtopq**2))**2)
     .       *(3d0*cosb**2*Ytq**2+(3d0*cosb**2+1d0)*Ybq**2))*vuq**4
     .    +3d0/512d0/Pi**4*Ybq**4*(dlog(QSTSB/mtopq**2)**2
     .      *(64d0*Pi*ALSQ-2d0/3d0*g1q+3d0*sinb**2*Ytq**2
     .       -3d0*cosb**2*Ybq**2)+
     .       (dlog(MA2/mtopq**2)**2-dlog(QSTSB/mtopq**2)**2)
     .       *(3d0*sinb**2*Ybq**2+(3d0*sinb**2+1d0)*Ytq**2))*vdq**4
     .    -(MCH2(1)**2*(dlog(MCH2(1)/QSTSB)-1.5d0)
     .     +MCH2(2)**2*(dlog(MCH2(2)/QSTSB)-1.5d0))/16d0/Pi**2
     .    -(MNEU(1)**2*(dlog(MNEU(1)/QSTSB)-1.5d0)
     .     +MNEU(2)**2*(dlog(MNEU(2)/QSTSB)-1.5d0)
     .     +MNEU(3)**2*(dlog(MNEU(3)/QSTSB)-1.5d0)
     .     +MNEU(4)**2*(dlog(MNEU(4)/QSTSB)-1.5d0)
     .     +MNEU(5)**2*(dlog(MNEU(5)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +MHC**2*(dlog(MHC/QSTSB)-1.5d0)/32d0/Pi**2
     .    +(MH0(1)**2*(dlog(MH0(1)/QSTSB)-1.5d0)
     .     +MH0(2)**2*(dlog(MH0(2)/QSTSB)-1.5d0)
     .     +MH0(3)**2*(dlog(MH0(3)/QSTSB)-1.5d0)
     .     +MH0(4)**2*(dlog(MH0(4)/QSTSB)-1.5d0)
     .     +MH0(5)**2*(dlog(MH0(5)/QSTSB)-1.5d0))/64d0/Pi**2

c      Z3-symmetric minimum (or with vanishing vev's for the Z3-violating case)

      Taux=MHuS+MHdS
      Daux=MHuS-MHdS
      aux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
      MHD1aux=(Taux-dsqrt(Daux**2+4d0*aux))/2d0
      MHD2aux=(Taux+dsqrt(Daux**2+4d0*aux))/2d0

      Taux=2d0*(MSS+MUPQ**2)
      Daux=MSPQ**2+(2d0*k*XIFQ)**2+4d0*k*XIFQ*MSPQ*DDCOS(phi02-phiF)
      MS1aux=(Taux-dsqrt(Daux))/2d0
      MS2aux=(Taux+dsqrt(Daux))/2d0

      VZ3=(3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSU3**2*(dlog(MSU3/QSTSB)-1.5d0)
     .     +3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSD3**2*(dlog(MSD3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +MSE3**2*(dlog(MSE3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSU1**2*(dlog(MSU1/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSD1**2*(dlog(MSD1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0)
     .     +2d0*MSE1**2*(dlog(MSE1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0))/32d0/Pi**2
     .    -(M2r**4*(dlog(M2r**2/QSTSB)-1.5d0))/16d0/Pi**2
     .    -(M1r**4*(dlog(M1r**2/QSTSB)-1.5d0)
     .     +M2r**4*(dlog(M2r**2/QSTSB)-1.5d0)
     .     +mupsi**4*(dlog(Max(mupsi**2,1d0)/QSTSB)-1.5d0))/32d0/Pi**2
     .     +(4d0*MHD1aux**2*(dlog(Max(MHD1aux,MZ2)/QSTSB)-1.5d0)
     .      +4d0*MHD2aux**2*(dlog(Max(MHD2aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS1aux**2*(dlog(Max(MS1aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS2aux**2*(dlog(Max(MS2aux,MZ2)/QSTSB)-1.5d0))
     .                                             /64d0/Pi**2

      IF(VZ3.lt.VEW)PROB(28)=DDIM(-1d-2,(VZ3-VEW)/Max(dabs(VEW),1d-10))

c      EW-symmetric minimum

      Taux=MHuS+MHdS+2d0*muq**2
      Daux=MHuS-MHdS
      aux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
     . +(MUPQ*muq)**2+2d0*M3HQ*MUPQ*muq*DDCOS(PhiP-Phi01-Phi3q)
     . +2d0*l*XIFQ*MUPQ*muq*DDCOS(PhiP-PhiF)
     . +muq**2*(Alcos1**2+IAL**2)+(k*muq**2/l)**2
     . +2d0*muq**3*k/l*(Alcos1*DDCOS(phi0)+IAL*DDSIN(phi0))
     . +2d0*muq*Alcos1*(M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-PhiF)
     .                  +muq*MUPQ*DDCOS(phi01-phiP))
     . +2d0*muq*IAL*(M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-PhiF)
     .                  +muq*MUPQ*DDSIN(phi01-phiP))
     . +2d0*k/l*muq**2*(M3HQ*DDCOS(phi3q-phi0)+l*XIFQ*DDCOS(phi02-PhiF)
     .                  +muq*MUPQ*DDCOS(phi02-phiP))
      MHD1aux=(Taux-dsqrt(Daux**2+4d0*aux))/2d0
      MHD2aux=(Taux+dsqrt(Daux**2+4d0*aux))/2d0

      Taux=2d0*(MSS+MUPQ**2)
     .     +4d0*((k*muq/l)**2+k/l*MUPQ*muq*DDCOS(phi02-phiP))
      Daux=2d0*MSPQ*DDCOS(phiSPq)+4d0*k*XIFQ*DDCOS(phi02-phiF)
     . +8d0*(k/l*muq)**2+4d0*k/l*muq*(Akcos2+2d0*MUPQ*DDCOS(phi02-phiP))
      aux=MSPQ*DDSIN(PhiSPq)+2d0*k*XIFQ*DDSIN(phi02-phiF)
     . +2d0*k/l*muq*(IAk+MUPQ*DDSIN(phi02-phiP))
      MS1aux=(Taux-dsqrt(Daux**2+4d0*aux**2))/2d0
      MS2aux=(Taux+dsqrt(Daux**2+4d0*aux**2))/2d0

      Vgauge=MSS*(muq/l)**2+2d0/3d0*k*Akcos2*(muq/l)**3+k**2*(muq/l)**4
     .    +2d0*(XISQ*DDCOS(phiSq)+XIFQ*MUPQ*DDCOS(phiP-phiF))*muq/l
     .    +2d0*(MSPQ/2d0*DDCOS(phiSPq)+k*XIFQ*DDCOS(phi02-phiF))
     .     *(muq/l)**2+2d0*k*MUPQ*DDCOS(phi02-phiP)*(muq/l)**3
     .    +(MUPQ*muq/l)**2
     .    +(3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSU3**2*(dlog(MSU3/QSTSB)-1.5d0)
     .     +3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSD3**2*(dlog(MSD3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +MSE3**2*(dlog(MSE3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSU1**2*(dlog(MSU1/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSD1**2*(dlog(MSD1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0)
     .     +2d0*MSE1**2*(dlog(MSE1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0))/32d0/Pi**2
     .    -(M2r**4*(dlog(M2r**2/QSTSB)-1.5d0)
     .     +mur**4*(dlog(mur**2/QSTSB)-1.5d0))/16d0/Pi**2
     .    -(M1r**4*(dlog(M1r**2/QSTSB)-1.5d0)
     .     +M2r**4*(dlog(M2r**2/QSTSB)-1.5d0)
     .     +2d0*mur**4*(dlog(mur**2/QSTSB)-1.5d0)
     .     +mupsi**4*(dlog(Max(mupsi**2,1d0)/QSTSB)-1.5d0))/32d0/Pi**2
     .     +(4d0*MHD1aux**2*(dlog(Max(MHD1aux,MZ2)/QSTSB)-1.5d0)
     .      +4d0*MHD2aux**2*(dlog(Max(MHD2aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS1aux**2*(dlog(Max(MS1aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS2aux**2*(dlog(Max(MS2aux,MZ2)/QSTSB)-1.5d0))
     .                                             /64d0/Pi**2

      Aroot=4d0*k**2
      Broot=2d0*k*Akcos2+6d0*k*MUPQ*DDCOS(phi02-phiP)
      Croot=2d0*MSS+4d0*(MSPQ/2d0*DDCOS(phiSPq)
     .     +k*XIFQ*DDCOS(phi02-phiF))+2d0*MUPQ**2
      Droot=2d0*(XISQ*DDCOS(phiSq)+XIFQ*MUPQ*DDCOS(phiP-phiF))
      root1=0d0
      root2=0d0
      root3=0d0

        IF(Aroot.ne.0d0)THEN
      Proot=-(Broot/Aroot)**2/3d0+Croot/Aroot
      Qroot=Broot/27d0/Aroot*(2d0*(Broot/Aroot)**2-9d0*Croot/Aroot)
     .     +Droot/Aroot
      Delt=-4d0*Proot**3-27d0*Qroot**2
         IF(Delt.gt.0d0)THEN
      root1=-Broot/3d0/Aroot+2d0*dsqrt(-Proot/3d0)
     .      *DDCOS(DACOS(-Qroot/2d0*dsqrt(27d0/(-Proot**3))))
      root2=-Broot/3d0/Aroot+2d0*dsqrt(-Proot/3d0)
     .      *DDCOS(DACOS(-Qroot/2d0*dsqrt(27d0/(-Proot**3)))+2d0/3d0*Pi)
      root3=-Broot/3d0/Aroot+2d0*dsqrt(-Proot/3d0)
     .      *DDCOS(DACOS(-Qroot/2d0*dsqrt(27d0/(-Proot**3)))+4d0/3d0*Pi)
         ELSEIF(Delt.eq.0d0)THEN
      root1=-Broot/3d0/Aroot
         ELSEIF(Delt.lt.0d0)THEN
      root1=-Broot/3d0/Aroot
     .     +((-Qroot+dsqrt(-Delt/27d0))/2d0)**(1d0/3d0)
     .     +((-Qroot-dsqrt(-Delt/27d0))/2d0)**(1d0/3d0)
         ENDIF
        ELSE
         IF(Broot.ne.0d0)THEN
      Proot=Croot/Broot
      Qroot=Droot/Broot
      Delt=Proot**2-4d0*Qroot
          IF(Delt.gt.0d0)THEN
      root1=(-Proot-dsqrt(Delt))/2d0
      root2=(-Proot+dsqrt(Delt))/2d0
          ELSEIF(Delt.eq.0d0)THEN
      root1=-Proot/2d0
          ENDIF
         ELSE
          IF(Croot.ne.0d0)root1=-Droot/Croot
         ENDIF
        ENDIF

        If(root1.ne.0d0)then
      Taux=MHuS+MHdS+2d0*(l*root1)**2
      Daux=MHuS-MHdS
      aux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
     . +(MUPQ*l*root1)**2+2d0*M3HQ*MUPQ*l*root1*DDCOS(PhiP-Phi01-Phi3q)
     . +2d0*l*XIFQ*MUPQ*l*root1*DDCOS(PhiP-PhiF)
     . +(l*root1)**2*(Alcos1**2+IAL**2)+(k*l*root1**2)**2
     . +2d0*k*l**2*root1**3*(Alcos1*DDCOS(phi0)+IAL*DDSIN(phi0))
     . +2d0*l*root1*Alcos1*(M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-PhiF)
     .                  +l*root1*MUPQ*DDCOS(phi01-phiP))
     . +2d0*l*root1*IAL*(M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-PhiF)
     .                  +l*root1*MUPQ*DDSIN(phi01-phiP))
     . +2d0*k*l*root1**2*(M3HQ*DDCOS(phi3q-phi0)
     .                  +l*XIFQ*DDCOS(phi02-PhiF)
     .                  +l*root1*MUPQ*DDCOS(phi02-phiP))
      MHD1aux=(Taux-dsqrt(Daux**2+4d0*aux))/2d0
      MHD2aux=(Taux+dsqrt(Daux**2+4d0*aux))/2d0

      Taux=2d0*(MSS+MUPQ**2)
     .     +4d0*((k*root1)**2+k*MUPQ*root1*DDCOS(phi02-phiP))
      Daux=2d0*MSPQ*DDCOS(phiSPq)+4d0*k*XIFQ*DDCOS(phi02-phiF)
     . +8d0*(k*root1)**2+4d0*k*root1*(Akcos2+2d0*MUPQ*DDCOS(phi02-phiP))
      aux=MSPQ*DDSIN(PhiSPq)+2d0*k*XIFQ*DDSIN(phi02-phiF)
     . +2d0*k*root1*(IAk+MUPQ*DDSIN(phi02-phiP))
      MS1aux=(Taux-dsqrt(Daux**2+4d0*aux**2))/2d0
      MS2aux=(Taux+dsqrt(Daux**2+4d0*aux**2))/2d0

      Vgauge=min(Vgauge,
     .     MSS*(root1)**2+2d0/3d0*k*Akcos2*(root1)**3+k**2*(root1)**4
     .    +2d0*(XISQ*DDCOS(phiSq)+XIFQ*MUPQ*DDCOS(phiP-phiF))*root1
     .    +2d0*(MSPQ/2d0*DDCOS(phiSPq)+k*XIFQ*DDCOS(phi02-phiF))
     .     *(root1)**2+2d0*k*MUPQ*DDCOS(phi02-phiP)*(root1)**3
     .    +(MUPQ*root1)**2
     .    +(3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSU3**2*(dlog(MSU3/QSTSB)-1.5d0)
     .     +3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSD3**2*(dlog(MSD3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +MSE3**2*(dlog(MSE3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSU1**2*(dlog(MSU1/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSD1**2*(dlog(MSD1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0)
     .     +2d0*MSE1**2*(dlog(MSE1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0))/32d0/Pi**2
     .    -(M2r**4*(dlog(M2r**2/QSTSB)-1.5d0)
     .     +(l*root1)**4*(dlog((l*root1)**2/QSTSB)-1.5d0))/16d0/Pi**2
     .    -(M1r**4*(dlog(M1r**2/QSTSB)-1.5d0)
     .     +M2r**4*(dlog(M2r**2/QSTSB)-1.5d0)
     .     +2d0*(l*root1)**4*(dlog((l*root1)**2/QSTSB)-1.5d0)
     .     +(mupsi+2d0*k/l*root1)**4
     .         *(dlog(Max((mupsi+2d0*k/l*root1)**2,1d0)/QSTSB)-1.5d0)
     .                                                   )/32d0/Pi**2
     .     +(4d0*MHD1aux**2*(dlog(Max(MHD1aux,MZ2)/QSTSB)-1.5d0)
     .      +4d0*MHD2aux**2*(dlog(Max(MHD2aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS1aux**2*(dlog(Max(MS1aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS2aux**2*(dlog(Max(MS2aux,MZ2)/QSTSB)-1.5d0))
     .                                             /64d0/Pi**2)
        Endif
        If(root2.ne.0d0)then
      Taux=MHuS+MHdS+2d0*(l*root2)**2
      Daux=MHuS-MHdS
      aux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
     . +(MUPQ*l*root2)**2+2d0*M3HQ*MUPQ*l*root2*DDCOS(PhiP-Phi01-Phi3q)
     . +2d0*l*XIFQ*MUPQ*l*root2*DDCOS(PhiP-PhiF)
     . +(l*root2)**2*(Alcos1**2+IAL**2)+(k*l*root2**2)**2
     . +2d0*k*l**2*root2**3*(Alcos1*DDCOS(phi0)+IAL*DDSIN(phi0))
     . +2d0*l*root2*Alcos1*(M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-PhiF)
     .                  +l*root2*MUPQ*DDCOS(phi01-phiP))
     . +2d0*l*root2*IAL*(M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-PhiF)
     .                  +l*root2*MUPQ*DDSIN(phi01-phiP))
     . +2d0*k*l*root2**2*(M3HQ*DDCOS(phi3q-phi0)
     .                  +l*XIFQ*DDCOS(phi02-PhiF)
     .                  +l*root2*MUPQ*DDCOS(phi02-phiP))
      MHD1aux=(Taux-dsqrt(Daux**2+4d0*aux))/2d0
      MHD2aux=(Taux+dsqrt(Daux**2+4d0*aux))/2d0

      Taux=2d0*(MSS+MUPQ**2)
     .     +4d0*((k*root2)**2+k*MUPQ*root2*DDCOS(phi02-phiP))
      Daux=2d0*MSPQ*DDCOS(phiSPq)+4d0*k*XIFQ*DDCOS(phi02-phiF)
     . +8d0*(k*root2)**2+4d0*k*root2*(Akcos2+2d0*MUPQ*DDCOS(phi02-phiP))
      aux=MSPQ*DDSIN(PhiSPq)+2d0*k*XIFQ*DDSIN(phi02-phiF)
     . +2d0*k*root2*(IAk+MUPQ*DDSIN(phi02-phiP))
      MS1aux=(Taux-dsqrt(Daux**2+4d0*aux**2))/2d0
      MS2aux=(Taux+dsqrt(Daux**2+4d0*aux**2))/2d0

      Vgauge=min(Vgauge,
     .     MSS*(root2)**2+2d0/3d0*k*Akcos2*(root2)**3+k**2*(root2)**4
     .    +2d0*(XISQ*DDCOS(phiSq)+XIFQ*MUPQ*DDCOS(phiP-phiF))*root2
     .    +2d0*(MSPQ/2d0*DDCOS(phiSPq)+k*XIFQ*DDCOS(phi02-phiF))
     .     *(root2)**2+2d0*k*MUPQ*DDCOS(phi02-phiP)*(root2)**3
     .    +(MUPQ*root2)**2
     .    +(3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSU3**2*(dlog(MSU3/QSTSB)-1.5d0)
     .     +3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSD3**2*(dlog(MSD3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +MSE3**2*(dlog(MSE3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSU1**2*(dlog(MSU1/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSD1**2*(dlog(MSD1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0)
     .     +2d0*MSE1**2*(dlog(MSE1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0))/32d0/Pi**2
     .    -(M2r**4*(dlog(M2r**2/QSTSB)-1.5d0)
     .     +(l*root2)**4*(dlog((l*root2)**2/QSTSB)-1.5d0))/16d0/Pi**2
     .    -(M1r**4*(dlog(M1r**2/QSTSB)-1.5d0)
     .     +M2r**4*(dlog(M2r**2/QSTSB)-1.5d0)
     .     +2d0*(l*root2)**4*(dlog((l*root2)**2/QSTSB)-1.5d0)
     .     +(mupsi+2d0*k/l*root2)**4
     .         *(dlog(Max((mupsi+2d0*k/l*root2)**2,1d0)/QSTSB)-1.5d0)
     .                                                   )/32d0/Pi**2
     .     +(4d0*MHD1aux**2*(dlog(Max(MHD1aux,MZ2)/QSTSB)-1.5d0)
     .      +4d0*MHD2aux**2*(dlog(Max(MHD2aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS1aux**2*(dlog(Max(MS1aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS2aux**2*(dlog(Max(MS2aux,MZ2)/QSTSB)-1.5d0))
     .                                             /64d0/Pi**2)
        Endif
        If(root3.ne.0d0)then
      Taux=MHuS+MHdS+2d0*(l*root3)**2
      Daux=MHuS-MHdS
      aux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
     . +(MUPQ*l*root3)**2+2d0*M3HQ*MUPQ*l*root3*DDCOS(PhiP-Phi01-Phi3q)
     . +2d0*l*XIFQ*MUPQ*l*root3*DDCOS(PhiP-PhiF)
     . +(l*root3)**2*(Alcos1**2+IAL**2)+(k*l*root3**2)**2
     . +2d0*k*l**2*root3**3*(Alcos1*DDCOS(phi0)+IAL*DDSIN(phi0))
     . +2d0*l*root3*Alcos1*(M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-PhiF)
     .                  +l*root3*MUPQ*DDCOS(phi01-phiP))
     . +2d0*l*root3*IAL*(M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-PhiF)
     .                  +l*root3*MUPQ*DDSIN(phi01-phiP))
     . +2d0*k*l*root3**2*(M3HQ*DDCOS(phi3q-phi0)
     .                  +l*XIFQ*DDCOS(phi02-PhiF)
     .                  +l*root3*MUPQ*DDCOS(phi02-phiP))
      MHD1aux=(Taux-dsqrt(Daux**2+4d0*aux))/2d0
      MHD2aux=(Taux+dsqrt(Daux**2+4d0*aux))/2d0

      Taux=2d0*(MSS+MUPQ**2)
     .     +4d0*((k*root3)**2+k*MUPQ*root3*DDCOS(phi02-phiP))
      Daux=2d0*MSPQ*DDCOS(phiSPq)+4d0*k*XIFQ*DDCOS(phi02-phiF)
     . +8d0*(k*root3)**2+4d0*k*root3*(Akcos2+2d0*MUPQ*DDCOS(phi02-phiP))
      aux=MSPQ*DDSIN(PhiSPq)+2d0*k*XIFQ*DDSIN(phi02-phiF)
     . +2d0*k*root3*(IAk+MUPQ*DDSIN(phi02-phiP))
      MS1aux=(Taux-dsqrt(Daux**2+4d0*aux**2))/2d0
      MS2aux=(Taux+dsqrt(Daux**2+4d0*aux**2))/2d0
      Vgauge=min(Vgauge,
     .     MSS*(root3)**2+2d0/3d0*k*Akcos2*(root3)**3+k**2*(root3)**4
     .    +2d0*(XISQ*DDCOS(phiSq)+XIFQ*MUPQ*DDCOS(phiP-phiF))*root3
     .    +2d0*(MSPQ/2d0*DDCOS(phiSPq)+k*XIFQ*DDCOS(phi02-phiF))
     .     *(root3)**2+2d0*k*MUPQ*DDCOS(phi02-phiP)*(root3)**3
     .    +(MUPQ*root3)**2
     .    +(3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSU3**2*(dlog(MSU3/QSTSB)-1.5d0)
     .     +3d0*MSQ3**2*(dlog(MSQ3/QSTSB)-1.5d0)
     .     +3d0*MSD3**2*(dlog(MSD3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +MSE3**2*(dlog(MSE3/QSTSB)-1.5d0)
     .     +MSL3**2*(dlog(MSL3/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSU1**2*(dlog(MSU1/QSTSB)-1.5d0)
     .     +6d0*MSQ1**2*(dlog(MSQ1/QSTSB)-1.5d0)
     .     +6d0*MSD1**2*(dlog(MSD1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0)
     .     +2d0*MSE1**2*(dlog(MSE1/QSTSB)-1.5d0)
     .     +2d0*MSL1**2*(dlog(MSL1/QSTSB)-1.5d0))/32d0/Pi**2
     .    -(M2r**4*(dlog(M2r**2/QSTSB)-1.5d0)
     .     +(l*root3)**4*(dlog((l*root3)**2/QSTSB)-1.5d0))/16d0/Pi**2
     .    -(M1r**4*(dlog(M1r**2/QSTSB)-1.5d0)
     .     +M2r**4*(dlog(M2r**2/QSTSB)-1.5d0)
     .     +2d0*(l*root3)**4*(dlog((l*root3)**2/QSTSB)-1.5d0)
     .     +(mupsi+2d0*k/l*root3)**4
     .         *(dlog(Max((mupsi+2d0*k/l*root3)**2,1d0)/QSTSB)-1.5d0)
     .                                                   )/32d0/Pi**2
     .     +(4d0*MHD1aux**2*(dlog(Max(MHD1aux,MZ2)/QSTSB)-1.5d0)
     .      +4d0*MHD2aux**2*(dlog(Max(MHD2aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS1aux**2*(dlog(Max(MS1aux,MZ2)/QSTSB)-1.5d0)
     .      +2d0*MS2aux**2*(dlog(Max(MS2aux,MZ2)/QSTSB)-1.5d0))
     .                                             /64d0/Pi**2)
        Endif

      IF(Vgauge.lt.VEW)PROB(28)=Max(PROB(28),
     .                   DDIM(-1d-2,(Vgauge-VEW)/Max(dabs(VEW),1d-10)))

c      Minimum with <S>=<Hd>=0

      Taux=MHuS+MHdS+g2q/2d0*vuq**2
      Daux=MHuS-MHdS+g1q/2d0*vuq**2
      aux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
      MHD1aux=(Taux-dsqrt(Daux**2+4d0*aux))/2d0
      MHD2aux=(Taux+dsqrt(Daux**2+4d0*aux))/2d0

      MH2(1,1)=MHuS+3d0/4d0*(g1q+g2q)*vuq**2
      MH2(1,2)=-M3HQ*DDCOS(phi3q)-l*XIFQ*DDCOS(phi01-phiF)
      MH2(2,1)=MH2(1,2)
      MH2(1,3)=0d0
      MH2(3,1)=MH2(1,3)
      MH2(2,2)=MHdS+(l**2-(g1q+g2q)/4d0)*vuq**2
      MH2(2,3)=-l*(Alcos1+MUPQ*DDCOS(phi01-phiP))*vuq
      MH2(3,2)=MH2(2,3)
      MH2(3,3)=MSS+MUPQ**2+(l*vuq)**2+MSPQ*DDCOS(phiSPq)
     .                               +2d0*k*XIFQ*DDCOS(phi02-phiF)
      MH2(1,4)=0d0
      MH2(4,1)=MH2(1,4)
      MH2(2,4)=M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
      MH2(4,2)=MH2(2,4)
      MH2(3,4)=0d0
      MH2(4,3)=MH2(3,4)
      MH2(4,4)=MHuS+(g1q+g2q)/4d0*vuq**2
      MH2(1,5)=M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
      MH2(5,1)=MH2(1,5)
      MH2(2,5)=0d0
      MH2(5,2)=MH2(2,5)
      MH2(3,5)=l*(IAl+MUPQ*DDSIN(phi01-phiP))*vuq
      MH2(5,3)=MH2(3,5)
      MH2(4,5)=M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-phiF)
      MH2(5,4)=MH2(4,5)
      MH2(5,5)=MHdS+(l**2-(g1q+g2q)/4d0)*vuq**2
      MH2(1,6)=0d0
      MH2(6,1)=MH2(1,6)
      MH2(2,6)=l*(IAl+MUPQ*DDSIN(phi01-phiP))*vuq
      MH2(6,2)=MH2(2,6)
      MH2(3,6)=-MSPQ*DDSIN(phiSPq)-2d0*k*XIFQ*DDSIN(phi02-phiF)
      MH2(6,3)=MH2(3,6)
      MH2(4,6)=0d0
      MH2(6,4)=MH2(4,6)
      MH2(5,6)=l*(Alcos1+MUPQ*DDCOS(phi01-phiP))*vuq
      MH2(6,5)=MH2(4,5)
      MH2(6,6)=MSS+MUPQ**2+l**2*vuq**2-MSPQ*DDCOS(phiSPq)
     .                               -2d0*k*XIFQ*DDCOS(phi02-phiF)

      CALL DIAGN(6,MH2,VALPH,VECPH,1d-10)

        MSTaux1=(MSQ3+MSU3+2d0*mtopq**2-(g1q+g2q)/4d0*vuq**2
     . -dsqrt((MSQ3-MSU3+(5d0/3d0*g1q-g2q)/4d0*vuq**2)**2
     .        +4d0*(Ytq*vuq*ATP)**2))/2d0
        MSTaux2=(MSQ3+MSU3+2d0*mtopq**2-(g1q+g2q)/4d0*vuq**2
     . +dsqrt((MSQ3-MSU3+(5d0/3d0*g1q-g2q)/4d0*vuq**2)**2
     .        +4d0*(Ytq*vuq*ATP)**2))/2d0
        MSBaux1=MSQ3+(g1q/3d0+g2q)/4d0*vuq**2
        MSBaux2=MSD3+g1q/6d0*vuq**2
        MSLaux1=MSL3+(-g1q+g2q)/4d0*vuq**2
        MSLaux2=MSE3+g1q/2d0*vuq**2
        MSNTaux=MSL3-(g1q+g2q)/4d0*vuq**2
        MSUaux1=MSQ1+(g1q/3d0-g2q)/4d0*vuq**2
        MSUaux2=MSU1-g1q/3d0*vuq**2
        MSDaux1=MSQ1+(g1q/3d0+g2q)/4d0*vuq**2
        MSDaux2=MSD1+g1q/6d0*vuq**2
        MSEaux1=MSL1+(-g1q+g2q)/4d0*vuq**2
        MSEaux2=MSE1+g1q/2d0*vuq**2
        MSNEaux=MSL1-(g1q+g2q)/4d0*vuq**2

      Vvu=MHuS*vuq**2+(G1Q+G2Q)/8d0*vuq**4
     .    -3d0*mtopq**4*(dlog(mtopq**2/QSTSB)-1.5d0)/16d0/Pi**2
     .    +3d0*(2d0*(g2q/2d0*vuq**2)**2
     .           *(dlog(g2q*vuq**2/2d0/QSTSB)-1.5d0)
     .         +((g1q+g2q)*vuq**2/2d0)**2
     .           *(dlog((g1q+g2q)*vuq**2/2d0/QSTSB)-1.5d0))/64d0/Pi**2
     .    +(3d0*MSTaux1**2*(dlog(Max(MSTaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSTaux2**2*(dlog(Max(MSTaux2,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux1**2*(dlog(Max(MSBaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux2**2*(dlog(Max(MSBaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux1**2*(dlog(Max(MSLaux1,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux2**2*(dlog(Max(MSLaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSNTaux**2*(dlog(Max(MSNTaux,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSUaux1**2*(dlog(Max(MSUaux1,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSUaux2**2*(dlog(Max(MSUaux2,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSDaux1**2*(dlog(Max(MSDaux1,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSDaux2**2*(dlog(Max(MSDaux2,MZ2)/QSTSB)-1.5d0)
     .     +2d0*MSEaux1**2*(dlog(Max(MSEaux1,MZ2)/QSTSB)-1.5d0)
     .     +2d0*MSEaux2**2*(dlog(Max(MSEaux2,MZ2)/QSTSB)-1.5d0)
     . +2d0*MSNEaux**2*(dlog(Max(MSNEaux,MZ2)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +3d0/512d0/Pi**4*Ytq**4*((dlog(QSTSB/mtopq**2))**2
     .      *(64d0*Pi*ALSQ+4d0/3d0*g1q-3d0*sinb**2*Ytq**2
     .       +3d0*cosb**2*Ybq**2)+
     .       ((dlog(MA2/mtopq**2))**2-(dlog(QSTSB/mtopq**2))**2)
     .       *(3d0*cosb**2*Ytq**2+(3d0*cosb**2+1d0)*Ybq**2))*vuq**4
     .    -((M2r**2+g2q*vuq**2)**2
     .              *(dlog((M2r**2+g2q*vuq**2)/QSTSB)-1.5d0))/16d0/Pi**2
     .    -((M1r**2+g1q*vuq**2)**2
     .                   *(dlog((M1r**2+g1q*vuq**2)/QSTSB)-1.5d0)
     .     +(M2r**2+g2q*vuq**2)**2
     .                   *(dlog((M2r**2+g2q*vuq**2)/QSTSB)-1.5d0)
     .     +(mupsi**2+2d0*(l*vuq)**2)**2*
     . (dlog(Max(mupsi**2+2d0*(l*vuq)**2,1d0)/QSTSB)-1.5d0))/32d0/Pi**2
     .     +(MHD1aux**2*(dlog(Max(MHD1aux,MW2)/QSTSB)-1.5d0)
     .     +MHD2aux**2*(dlog(Max(MHD2aux,MW2)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +(VALPH(1)**2*(dlog(Max(MZ2,VALPH(1))/QSTSB)-1.5d0)
     .     +VALPH(2)**2*(dlog(Max(MZ2,VALPH(2))/QSTSB)-1.5d0)
     .     +VALPH(3)**2*(dlog(Max(MZ2,VALPH(3))/QSTSB)-1.5d0)
     .     +VALPH(4)**2*(dlog(Max(MZ2,VALPH(4))/QSTSB)-1.5d0)
     .     +VALPH(5)**2*(dlog(Max(MZ2,VALPH(5))/QSTSB)-1.5d0)
     .   +VALPH(6)**2*(dlog(Max(MZ2,VALPH(6))/QSTSB)-1.5d0))/64d0/Pi**2

        If(MHuS.lt.0d0)then
      aux=dsqrt(-4d0*MHuS/(G1Q+G2Q))

      Taux=MHuS-MHdS+g1q/2d0*aux**2
      Daux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
      MHD1aux=(MHuS+MHdS+g2q/2d0*aux**2-dsqrt(Taux**2+4d0*Daux))/2d0
      MHD2aux=(MHuS+MHdS+g2q/2d0*aux**2+dsqrt(Taux**2+4d0*Daux))/2d0

      MH2(1,1)=MHuS+3d0/4d0*(g1q+g2q)*aux**2
      MH2(1,2)=-M3HQ*DDCOS(phi3q)-l*XIFQ*DDCOS(phi01-phiF)
      MH2(2,1)=MH2(1,2)
      MH2(1,3)=0d0
      MH2(3,1)=MH2(1,3)
      MH2(2,2)=MHdS+(l**2-(g1q+g2q)/4d0)*aux**2
      MH2(2,3)=-l*(Alcos1+MUPQ*DDCOS(phi01-phiP))*aux
      MH2(3,2)=MH2(2,3)
      MH2(3,3)=MSS+MUPQ**2+(l*aux)**2+MSPQ*DDCOS(phiSPq)
     .                               +2d0*k*XIFQ*DDCOS(phi02-phiF)
      MH2(1,4)=0d0
      MH2(4,1)=MH2(1,4)
      MH2(2,4)=M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
      MH2(4,2)=MH2(2,4)
      MH2(3,4)=0d0
      MH2(4,3)=MH2(3,4)
      MH2(4,4)=MHuS+(g1q+g2q)/4d0*aux**2
      MH2(1,5)=M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
      MH2(5,1)=MH2(1,5)
      MH2(2,5)=0d0
      MH2(5,2)=MH2(2,5)
      MH2(3,5)=l*(IAl+MUPQ*DDSIN(phi01-phiP))*aux
      MH2(5,3)=MH2(3,5)
      MH2(4,5)=M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-phiF)
      MH2(5,4)=MH2(4,5)
      MH2(5,5)=MHdS+(l**2-(g1q+g2q)/4d0)*aux**2
      MH2(1,6)=0d0
      MH2(6,1)=MH2(1,6)
      MH2(2,6)=l*(IAl+MUPQ*DDSIN(phi01-phiP))*aux
      MH2(6,2)=MH2(2,6)
      MH2(3,6)=-MSPQ*DDSIN(phiSPq)-2d0*k*XIFQ*DDSIN(phi02-phiF)
      MH2(6,3)=MH2(3,6)
      MH2(4,6)=0d0
      MH2(6,4)=MH2(4,6)
      MH2(5,6)=l*(Alcos1+MUPQ*DDCOS(phi01-phiP))*aux
      MH2(6,5)=MH2(4,5)
      MH2(6,6)=MSS+MUPQ**2+l**2*aux**2-MSPQ*DDCOS(phiSPq)
     .                               -2d0*k*XIFQ*DDCOS(phi02-phiF)

      CALL DIAGN(6,MH2,VALPH,VECPH,1d-10)

      MSTaux1=(MSQ3+MSU3+2d0*(Ytq*aux)**2-(g1q+g2q)/4d0*aux**2
     . -dsqrt((MSQ3-MSU3+(5d0/3d0*g1q-g2q)/4d0*aux**2)**2
     .        +4d0*(Ytq*aux*ATP)**2))/2d0
      MSTaux2=(MSQ3+MSU3+2d0*(Ytq*aux)**2-(g1q+g2q)/4d0*aux**2
     . +dsqrt((MSQ3-MSU3+(5d0/3d0*g1q-g2q)/4d0*aux**2)**2
     .        +4d0*(Ytq*aux*ATP)**2))/2d0
      MSBaux1=MSQ3+(g1q/3d0+g2q)/4d0*aux**2
      MSBaux2=MSD3+g1q/6d0*aux**2
      MSLaux1=MSL3+(-g1q+g2q)/4d0*aux**2
      MSLaux2=MSE3+g1q/2d0*aux**2
      MSNTaux=MSL3-(g1q+g2q)/4d0*aux**2
      MSNTaux=MSL3-(g1q+g2q)/4d0*aux**2
      MSUaux1=MSQ1+(g1q/3d0-g2q)/4d0*aux**2
      MSUaux2=MSU1-g1q/3d0*aux**2
      MSDaux1=MSQ1+(g1q/3d0+g2q)/4d0*aux**2
      MSDaux2=MSD1+g1q/6d0*aux**2
      MSEaux1=MSL1+(-g1q+g2q)/4d0*aux**2
      MSEaux2=MSE1+g1q/2d0*aux**2
      MSNEaux=MSL1-(g1q+g2q)/4d0*aux**2

      Vvu=min(Vvu,MHuS*aux**2+(G1Q+G2Q)/8d0*aux**4
     .    -3d0*(Ytq*aux)**4*(dlog((Ytq*aux)**2/QSTSB)-1.5d0)/16d0/Pi**2
     .    +3d0*(2d0*(g2q/2d0*aux**2)**2
     .           *(dlog(g2q*aux**2/2d0/QSTSB)-1.5d0)
     .         +((g1q+g2q)*aux**2/2d0)**2
     .           *(dlog((g1q+g2q)*aux**2/2d0/QSTSB)-1.5d0))/64d0/Pi**2
     .    +(3d0*MSTaux1**2*(dlog(Max(MSTaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSTaux2**2*(dlog(Max(MSTaux2,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux1**2*(dlog(Max(MSBaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux2**2*(dlog(Max(MSBaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux1**2*(dlog(Max(MSLaux1,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux2**2*(dlog(Max(MSLaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSNTaux**2*(dlog(Max(MSNTaux,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSUaux1**2*(dlog(Max(MSUaux1,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSUaux2**2*(dlog(Max(MSUaux2,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSDaux1**2*(dlog(Max(MSDaux1,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSDaux2**2*(dlog(Max(MSDaux2,MZ2)/QSTSB)-1.5d0)
     .     +2d0*MSEaux1**2*(dlog(Max(MSEaux1,MZ2)/QSTSB)-1.5d0)
     .     +2d0*MSEaux2**2*(dlog(Max(MSEaux2,MZ2)/QSTSB)-1.5d0)
     . +2d0*MSNEaux**2*(dlog(Max(MSNEaux,MZ2)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +3d0/512d0/Pi**4*Ytq**4*((dlog(QSTSB/mtopq**2))**2
     .      *(64d0*Pi*ALSQ+4d0/3d0*g1q-3d0*sinb**2*Ytq**2
     .       +3d0*cosb**2*Ybq**2)+
     .       ((dlog(MA2/mtopq**2))**2-(dlog(QSTSB/mtopq**2))**2)
     .       *(3d0*cosb**2*Ytq**2+(3d0*cosb**2+1d0)*Ybq**2))*aux**4
     .    -((M2r**2+g2q*aux**2)**2
     .            *(dlog((M2r**2+g2q*aux**2)/QSTSB)-1.5d0))/16d0/Pi**2
     .    -((M1r**2+g1q*aux**2)**2
     .                   *(dlog((M1r**2+g1q*aux**2)/QSTSB)-1.5d0)
     .     +(M2r**2+g2q*aux**2)**2
     .                   *(dlog((M2r**2+g2q*aux**2)/QSTSB)-1.5d0)
     .     +(mupsi**2+2d0*(l*aux)**2)**2*
     . (dlog(Max(mupsi**2+2d0*(l*aux)**2,1d0)/QSTSB)-1.5d0))/32d0/Pi**2
     .     +(MHD1aux**2*(dlog(Max(MHD1aux,MW2)/QSTSB)-1.5d0)
     .     +MHD2aux**2*(dlog(Max(MHD2aux,MW2)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +(VALPH(1)**2*(dlog(Max(MZ2,VALPH(1))/QSTSB)-1.5d0)
     .     +VALPH(2)**2*(dlog(Max(MZ2,VALPH(2))/QSTSB)-1.5d0)
     .     +VALPH(3)**2*(dlog(Max(MZ2,VALPH(3))/QSTSB)-1.5d0)
     .     +VALPH(4)**2*(dlog(Max(MZ2,VALPH(4))/QSTSB)-1.5d0)
     .     +VALPH(5)**2*(dlog(Max(MZ2,VALPH(5))/QSTSB)-1.5d0)
     .   +VALPH(6)**2*(dlog(Max(MZ2,VALPH(6))/QSTSB)-1.5d0))/64d0/Pi**2)
        endif

      IF(Vvu.lt.VEW)PROB(28)=Max(PROB(28),
     .                      DDIM(-1d-2,(VEW-Vvu)/Max(dabs(VEW),1d-10)))

c      Minimum with <S>=<Hu>=0

      Taux=MHuS+MHdS+g2q/2d0*vdq**2
      Daux=MHuS-MHdS-g1q/2d0*vdq**2
      aux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
      MHD1aux=(Taux-dsqrt(Daux**2+4d0*aux))/2d0
      MHD2aux=(Taux+dsqrt(Daux**2+4d0*aux))/2d0

      MH2(1,1)=MHuS+(l**2-(g1q+g2q)/4d0)*vdq**2
      MH2(1,2)=-M3HQ*DDCOS(phi3q)-l*XIFQ*DDCOS(phi01-phiF)
      MH2(2,1)=MH2(1,2)
      MH2(1,3)=-l*(Alcos1+MUPQ*DDCOS(phi01-phiP))*vdq
      MH2(3,1)=MH2(1,3)
      MH2(2,2)=MHdS+3d0*(g1q+g2q)/4d0*vdq**2
      MH2(2,3)=0d0
      MH2(3,2)=MH2(2,3)
      MH2(3,3)=MSS+MUPQ**2+(l*vdq)**2+MSPQ*DDCOS(phiSPq)
     .                               +2d0*k*XIFQ*DDCOS(phi02-phiF)
      MH2(1,4)=0d0
      MH2(4,1)=MH2(1,4)
      MH2(2,4)=M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
      MH2(4,2)=MH2(2,4)
      MH2(3,4)=l*(IAl+MUPQ*DDSIN(phi01-phiP))*vdq
      MH2(4,3)=MH2(3,4)
      MH2(4,4)=MHuS+(l**2-(g1q+g2q)/4d0)*vdq**2
      MH2(1,5)=M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
      MH2(5,1)=MH2(1,5)
      MH2(2,5)=0d0
      MH2(5,2)=MH2(2,5)
      MH2(3,5)=0d0
      MH2(5,3)=MH2(3,5)
      MH2(4,5)=M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-phiF)
      MH2(5,4)=MH2(4,5)
      MH2(5,5)=MHdS+(g1q+g2q)/4d0*vdq**2
      MH2(1,6)=l*(IAl+MUPQ*DDSIN(phi01-phiP))*vdq
      MH2(6,1)=MH2(1,6)
      MH2(2,6)=0d0
      MH2(6,2)=MH2(2,6)
      MH2(3,6)=-MSPQ*DDSIN(phiSPq)-2d0*k*XIFQ*DDSIN(phi02-phiF)
      MH2(6,3)=MH2(3,6)
      MH2(4,6)=l*(Alcos1+MUPQ*DDCOS(phi01-phiP))*vdq
      MH2(6,4)=MH2(4,6)
      MH2(5,6)=0d0
      MH2(6,5)=MH2(4,5)
      MH2(6,6)=MSS+MUPQ**2+l**2*vdq**2-MSPQ*DDCOS(phiSPq)
     .                               -2d0*k*XIFQ*DDCOS(phi02-phiF)

      CALL DIAGN(6,MH2,VALPH,VECPH,1d-10)

        MSTaux1=MSQ3-(g1q/3d0-g2q)/4d0*vdq**2
        MSTaux2=MSU3+g1q/3d0*vdq**2
        MSBaux1=(MSQ3+MSD3+2d0*mbotq**2-(g1q+g2q)/4d0*vdq**2
     . -dsqrt((MSQ3-MSD3+(-g1q/3d0+g2q)/4d0*vdq**2)**2
     .        +4d0*(Ybq*vdq*ABP)**2))/2d0
        MSBaux2=(MSQ3+MSD3+2d0*mbotq**2-(g1q+g2q)/4d0*vdq**2
     . +dsqrt((MSQ3-MSD3+(-g1q/3d0+g2q)/4d0*vdq**2)**2
     .        +4d0*(Ybq*vdq*ABP)**2))/2d0
        MSLaux1=(MSL3+MSE3+2d0*(mtau*vdq/vd)**2-(g1q+g2q)/4d0*vdq**2
     . -dsqrt((MSL3-MSE3-(-3d0/2d0*g1q+g2q)/4d0*vdq**2)**2
     .        +4d0*(mtau/vd*vdq*ATAU)**2))/2d0
        MSLaux2=(MSL3+MSE3+2d0*(mtau*vdq/vd)**2-(g1q+g2q)/4d0*vdq**2
     . +dsqrt((MSL3-MSE3-(-3d0/2d0*g1q+g2q)/4d0*vdq**2)**2
     .        +4d0*(mtau/vd*vdq*ATAU)**2))/2d0
        MSNTaux=MSL3+(g1q+g2q)/4d0*vdq**2
        MSUaux1=MSQ1-(g1q/3d0-g2q)/4d0*vdq**2
        MSUaux2=MSU1+g1q/3d0*vdq**2
        MSDaux1=MSQ1-(g1q/3d0+g2q)/4d0*vdq**2
        MSDaux2=MSD1-g1q/6d0*vdq**2
        MSEaux1=MSL1-(-g1q+g2q)/4d0*vdq**2
        MSEaux2=MSE1-g1q/2d0*vdq**2
        MSNEaux=MSL1+(g1q+g2q)/4d0*vdq**2

      Vvd=MHdS*vdq**2+(G1Q+G2Q)/8d0*vdq**4
     .    -(3d0*mbotq**4*(dlog(mbotq**2/QSTSB)-1.5d0)
     .     +(mtau*vdq/vd)**4*(dlog(mtau**2/QSTSB)-1.5d0))/16d0/Pi**2
     .    +3d0*(2d0*(g2q*vdq**2/2d0)**2
     .           *(dlog(g2q*vdq**2/2d0/QSTSB)-1.5d0)
     .         +((g1q+g2q)/2d0*vdq**2)**2
     .           *(dlog((g1q+g2q)/2d0*vdq**2/QSTSB)-1.5d0))/64d0/Pi**2
     .    +(3d0*MSTaux1**2*(dlog(Max(MSTaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSTaux2**2*(dlog(Max(MSTaux2,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux1**2*(dlog(Max(MSBaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux2**2*(dlog(Max(MSBaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux1**2*(dlog(Max(MSLaux1,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux2**2*(dlog(Max(MSLaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSNTaux**2*(dlog(Max(MSNTaux,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSUaux1**2*(dlog(Max(MSUaux1,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSUaux2**2*(dlog(Max(MSUaux2,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSDaux1**2*(dlog(Max(MSDaux1,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSDaux2**2*(dlog(Max(MSDaux2,MZ2)/QSTSB)-1.5d0)
     .     +2d0*MSEaux1**2*(dlog(Max(MSEaux1,MZ2)/QSTSB)-1.5d0)
     .     +2d0*MSEaux2**2*(dlog(Max(MSEaux2,MZ2)/QSTSB)-1.5d0)
     . +2d0*MSNEaux**2*(dlog(Max(MSNEaux,MZ2)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +3d0/512d0/Pi**4*Ybq**4*(dlog(QSTSB/mtopq**2)**2
     .      *(64d0*Pi*ALSQ-2d0/3d0*g1q+3d0*sinb**2*Ytq**2
     .       -3d0*cosb**2*Ybq**2)+
     .       (dlog(MA2/mtopq**2)**2-dlog(QSTSB/mtopq**2)**2)
     .       *(3d0*sinb**2*Ybq**2+(3d0*sinb**2+1d0)*Ytq**2))*vdq**4
     .    -((M2r**2+g2q*vdq**2)**2
     .             *(dlog((M2r**2+g2q*vdq**2)/QSTSB)-1.5d0))/16d0/Pi**2
     .    -((M1r**2+g1q*vdq**2)**2
     .                   *(dlog((M1r**2+g1q*vdq**2)/QSTSB)-1.5d0)
     .     +(M2r**2+g2q*vdq**2)**2
     .                   *(dlog((M2r**2+g2q*vdq**2)/QSTSB)-1.5d0)
     .     +(mupsi**2+2d0*(l*vdq)**2)**2*
     . (dlog(Max(mupsi**2+2d0*(l*vuq)**2,1d0)/QSTSB)-1.5d0))/32d0/Pi**2
     .     +(MHD1aux**2*(dlog(Max(MHD1aux,MW2)/QSTSB)-1.5d0)
     .     +MHD2aux**2*(dlog(Max(MHD2aux,MW2)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +(VALPH(1)**2*(dlog(Max(MZ2,VALPH(1))/QSTSB)-1.5d0)
     .     +VALPH(2)**2*(dlog(Max(MZ2,VALPH(2))/QSTSB)-1.5d0)
     .     +VALPH(3)**2*(dlog(Max(MZ2,VALPH(3))/QSTSB)-1.5d0)
     .     +VALPH(4)**2*(dlog(Max(MZ2,VALPH(4))/QSTSB)-1.5d0)
     .     +VALPH(5)**2*(dlog(Max(MZ2,VALPH(5))/QSTSB)-1.5d0)
     .   +VALPH(6)**2*(dlog(Max(MZ2,VALPH(6))/QSTSB)-1.5d0))/64d0/Pi**2

        If(MHdS.lt.0d0)then
      aux=dsqrt(-4d0*MHdS/(G1Q+G2Q))

      Taux=MHuS-MHdS-g1q/2d0*aux**2
      Daux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
      MHD1aux=(MHuS+MHdS+g2q/2d0*aux**2-dsqrt(Taux**2+4d0*Daux))/2d0
      MHD2aux=(MHuS+MHdS+g2q/2d0*aux**2+dsqrt(Taux**2+4d0*Daux))/2d0

      MH2(1,1)=MHuS+(l**2-(g1q+g2q)/4d0)*aux**2
      MH2(1,2)=-M3HQ*DDCOS(phi3q)-l*XIFQ*DDCOS(phi01-phiF)
      MH2(2,1)=MH2(1,2)
      MH2(1,3)=-l*(Alcos1+MUPQ*DDCOS(phi01-phiP))*aux
      MH2(3,1)=MH2(1,3)
      MH2(2,2)=MHdS+3d0*(g1q+g2q)/4d0*aux**2
      MH2(2,3)=0d0
      MH2(3,2)=MH2(2,3)
      MH2(3,3)=MSS+MUPQ**2+(l*aux)**2+MSPQ*DDCOS(phiSPq)
     .                               +2d0*k*XIFQ*DDCOS(phi02-phiF)
      MH2(1,4)=0d0
      MH2(4,1)=MH2(1,4)
      MH2(2,4)=M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
      MH2(4,2)=MH2(2,4)
      MH2(3,4)=l*(IAl+MUPQ*DDSIN(phi01-phiP))*aux
      MH2(4,3)=MH2(3,4)
      MH2(4,4)=MHuS+(l**2-(g1q+g2q)/4d0)*aux**2
      MH2(1,5)=M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
      MH2(5,1)=MH2(1,5)
      MH2(2,5)=0d0
      MH2(5,2)=MH2(2,5)
      MH2(3,5)=0d0
      MH2(5,3)=MH2(3,5)
      MH2(4,5)=M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-phiF)
      MH2(5,4)=MH2(4,5)
      MH2(5,5)=MHdS+(g1q+g2q)/4d0*aux**2
      MH2(1,6)=l*(IAl+MUPQ*DDSIN(phi01-phiP))*aux
      MH2(6,1)=MH2(1,6)
      MH2(2,6)=0d0
      MH2(6,2)=MH2(2,6)
      MH2(3,6)=-MSPQ*DDSIN(phiSPq)-2d0*k*XIFQ*DDSIN(phi02-phiF)
      MH2(6,3)=MH2(3,6)
      MH2(4,6)=l*(Alcos1+MUPQ*DDCOS(phi01-phiP))*aux
      MH2(6,4)=MH2(4,6)
      MH2(5,6)=0d0
      MH2(6,5)=MH2(4,5)
      MH2(6,6)=MSS+MUPQ**2+l**2*aux**2-MSPQ*DDCOS(phiSPq)
     .                               -2d0*k*XIFQ*DDCOS(phi02-phiF)

      CALL DIAGN(6,MH2,VALPH,VECPH,1d-10)

      MSTaux1=MSQ3-(g1q/3d0-g2q)/4d0*aux**2
        MSTaux2=MSU3+g1q/3d0*aux**2
        MSBaux1=(MSQ3+MSD3+2d0*(Ybq*aux)**2-(g1q+g2q)/4d0*aux**2
     . -dsqrt((MSQ3-MSD3+(-g1q/3d0+g2q)/4d0*aux**2)**2
     .        +4d0*(Ybq*aux*ABP)**2))/2d0
        MSBaux2=(MSQ3+MSD3+2d0*(Ybq*aux)**2-(g1q+g2q)/4d0*aux**2
     . +dsqrt((MSQ3-MSD3+(-g1q/3d0+g2q)/4d0*aux**2)**2
     .        +4d0*(Ybq*aux*ABP)**2))/2d0
        MSLaux1=(MSL3+MSE3+2d0*(mtau*aux/vd)**2-(g1q+g2q)/4d0*aux**2
     . -dsqrt((MSL3-MSE3-(-3d0/2d0*g1q+g2q)/4d0*aux**2)**2
     .        +4d0*(mtau/vd*aux*ATAU)**2))/2d0
        MSLaux2=(MSL3+MSE3+2d0*(mtau*aux/vd)**2-(g1q+g2q)/4d0*aux**2
     . +dsqrt((MSL3-MSE3-(-3d0/2d0*g1q+g2q)/4d0*aux**2)**2
     .        +4d0*(mtau/vd*aux*ATAU)**2))/2d0
        MSNTaux=MSL3+(g1q+g2q)/4d0*aux**2
        MSUaux1=MSQ1-(g1q/3d0-g2q)/4d0*aux**2
        MSUaux2=MSU1+g1q/3d0*aux**2
        MSDaux1=MSQ1-(g1q/3d0+g2q)/4d0*aux**2
        MSDaux2=MSD1-g1q/6d0*aux**2
        MSEaux1=MSL1-(-g1q+g2q)/4d0*aux**2
        MSEaux2=MSE1-g1q/2d0*aux**2
        MSNEaux=MSL1+(g1q+g2q)/4d0*aux**2

      Vvd=min(Vvd,MHdS*aux**2+(G1Q+G2Q)/8d0*aux**4
     .    -(3d0*(Ybq*aux)**4*(dlog((Ybq*aux)**2/QSTSB)-1.5d0)
     .     +(mtau*aux/vd)**4*(dlog((mtau*aux/vd)**2/QSTSB)-1.5d0)
     .                                                    )/16d0/Pi**2
     .    +3d0*(2d0*(g2q*aux**2/2d0)**2
     .           *(dlog(g2q*aux**2/2d0/QSTSB)-1.5d0)
     .         +((g1q+g2q)/2d0*aux**2)**2
     .           *(dlog((g1q+g2q)/2d0*aux**2/QSTSB)-1.5d0))/64d0/Pi**2
     .    +(3d0*MSTaux1**2*(dlog(Max(MSTaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSTaux2**2*(dlog(Max(MSTaux2,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux1**2*(dlog(Max(MSBaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux2**2*(dlog(Max(MSBaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux1**2*(dlog(Max(MSLaux1,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux2**2*(dlog(Max(MSLaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSNTaux**2*(dlog(Max(MSNTaux,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSUaux1**2*(dlog(Max(MSUaux1,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSUaux2**2*(dlog(Max(MSUaux2,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSDaux1**2*(dlog(Max(MSDaux1,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSDaux2**2*(dlog(Max(MSDaux2,MZ2)/QSTSB)-1.5d0)
     .     +2d0*MSEaux1**2*(dlog(Max(MSEaux1,MZ2)/QSTSB)-1.5d0)
     .     +2d0*MSEaux2**2*(dlog(Max(MSEaux2,MZ2)/QSTSB)-1.5d0)
     . +2d0*MSNEaux**2*(dlog(Max(MSNEaux,MZ2)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +3d0/512d0/Pi**4*Ybq**4*(dlog(QSTSB/mtopq**2)**2
     .      *(64d0*Pi*ALSQ-2d0/3d0*g1q+3d0*sinb**2*Ytq**2
     .       -3d0*cosb**2*Ybq**2)+
     .       (dlog(MA2/mtopq**2)**2-dlog(QSTSB/mtopq**2)**2)
     .       *(3d0*sinb**2*Ybq**2+(3d0*sinb**2+1d0)*Ytq**2))*aux**4
     .    -((M2r**2+g2q*aux**2)**2
     .             *(dlog((M2r**2+g2q*aux**2)/QSTSB)-1.5d0))/16d0/Pi**2
     .    -((M1r**2+g1q*aux**2)**2
     .                   *(dlog((M1r**2+g1q*aux**2)/QSTSB)-1.5d0)
     .     +(M2r**2+g2q*aux**2)**2
     .                   *(dlog((M2r**2+g2q*aux**2)/QSTSB)-1.5d0)
     .     +(mupsi**2+2d0*(l*aux)**2)**2*
     . (dlog(Max(mupsi**2+2d0*(l*vuq)**2,1d0)/QSTSB)-1.5d0))/32d0/Pi**2
     .     +(MHD1aux**2*(dlog(Max(MHD1aux,MW2)/QSTSB)-1.5d0)
     .     +MHD2aux**2*(dlog(Max(MHD2aux,MW2)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +(VALPH(1)**2*(dlog(Max(MZ2,VALPH(1))/QSTSB)-1.5d0)
     .     +VALPH(2)**2*(dlog(Max(MZ2,VALPH(2))/QSTSB)-1.5d0)
     .     +VALPH(3)**2*(dlog(Max(MZ2,VALPH(3))/QSTSB)-1.5d0)
     .     +VALPH(4)**2*(dlog(Max(MZ2,VALPH(4))/QSTSB)-1.5d0)
     .     +VALPH(5)**2*(dlog(Max(MZ2,VALPH(5))/QSTSB)-1.5d0)
     .   +VALPH(6)**2*(dlog(Max(MZ2,VALPH(6))/QSTSB)-1.5d0))/64d0/Pi**2)
        endif

      IF(Vvd.lt.VEW)PROB(28)=Max(PROB(28),
     .                      DDIM(-1d-2,(Vvd-VEW)/Max(dabs(VEW),1d-10)))

c      Phase variations

      DO I=1,Nphi
      DO J=1,Nphi

      phiS=2d0*Pi*I/DFLOAT(Nphi+1)
      phiH=2d0*Pi*J/DFLOAT(Nphi+1)

        MSTaux1=(MSQ3+MSU3+2d0*mtopq**2-(g1q+g2q)/4d0*(vuq**2-vdq**2)
     . -dsqrt((MSQ3-MSU3+(5d0/3d0*g1q-g2q)/4d0*(vuq**2-vdq**2))**2
     .        +4d0*Ytq**2*(ATP**2*vuq**2+muq**2*vdq**2
     .              -2d0*ATP*muq*vuq*vdq*DDCOS(phiAT+phiS+phiH))))/2d0
        MSTaux2=(MSQ3+MSU3+2d0*mtopq**2-(g1q+g2q)/4d0*(vuq**2-vdq**2)
     . +dsqrt((MSQ3-MSU3+(5d0/3d0*g1q-g2q)/4d0*(vuq**2-vdq**2))**2
     .        +4d0*Ytq**2*(ATP**2*vuq**2+muq**2*vdq**2
     .              -2d0*ATP*muq*vuq*vdq*DDCOS(phiAT+phiS+phiH))))/2d0
        MSBaux1=(MSQ3+MSD3+2d0*mbotq**2-(g1q+g2q)/4d0*(vdq**2-vuq**2)
     . -dsqrt((MSQ3-MSD3+(-g1q/3d0+g2q)/4d0*(vdq**2-vuq**2))**2
     .        +4d0*Ybq**2*(ABP**2*vdq**2+muq**2*vuq**2
     .              -2d0*ABP*muq*vuq*vdq*DDCOS(phiAB+phiS+phiH))))/2d0
        MSBaux2=(MSQ3+MSD3+2d0*mbotq**2-(g1q+g2q)/4d0*(vdq**2-vuq**2)
     . +dsqrt((MSQ3-MSD3+(-g1q/3d0+g2q)/4d0*(vdq**2-vuq**2))**2
     .        +4d0*Ybq**2*(ABP**2*vdq**2+muq**2*vuq**2
     .              -2d0*ABP*muq*vuq*vdq*DDCOS(phiAB+phiS+phiH))))/2d0
        MSLaux1=(MSL3+MSE3+2d0*(mtau*vdq/vd)**2
     .                                 -(g1q+g2q)/4d0*(vdq**2-vuq**2)
     . -dsqrt((MSL3-MSE3-(-3d0/2d0*g1q+g2q)/4d0*(vdq**2-vuq**2))**2
     .        +4d0*(mtau/vd)**2*(ATAU**2*vdq**2+muq**2*vuq**2
     .           -2d0*ATAU*muq*vuq*vdq*DDCOS(PhiATAU+PhiS+PhiH))))/2d0
        MSLaux2=(MSL3+MSE3+2d0*(mtau*vdq/vd)**2
     .                                 -(g1q+g2q)/4d0*(vdq**2-vuq**2)
     . +dsqrt((MSL3-MSE3-(-3d0/2d0*g1q+g2q)/4d0*(vdq**2-vuq**2))**2
     .        +4d0*(mtau/vd)**2*(ATAU**2*vdq**2+muq**2*vuq**2
     .           -2d0*ATAU*muq*vuq*vdq*DDCOS(PhiATAU+PhiS+PhiH))))/2d0
        MSNTaux=MSNT2

      Taux=M2r**2+mur**2+g2q*(vuq**2+vdq**2)
      Daux=M2r**2-mur**2+g2q*(vuq**2-vdq**2)
      aux=g2q*((M2r*vuq)**2+(mur*vdq)**2
     .         +2d0*mur*M2r*vuq*vdq*DDCOS(phiM2-phi01-phiS-phiH))
        MCH1aux=(Taux-dsqrt(Daux**2+4d0*aux))/2d0
        MCH2aux=(Taux+dsqrt(Daux**2+4d0*aux))/2d0

      MN2(1,1)=M1r**2+g1q/2d0*(vuq**2+vdq**2)
      MN2(6,6)=M1r**2+g1q/2d0*(vuq**2+vdq**2)
      MN2(6,1)=0.d0
      MN2(1,6)=0.d0
      MN2(1,2)=-dsqrt(g1q*g2q)/2d0*(vuq**2+vdq**2)
      MN2(2,1)=-dsqrt(g1q*g2q)/2d0*(vuq**2+vdq**2)
      MN2(6,7)=-dsqrt(g1q*g2q)/2d0*(vuq**2+vdq**2)
      MN2(7,6)=-dsqrt(g1q*g2q)/2d0*(vuq**2+vdq**2)
      MN2(1,7)=0.d0
      MN2(2,6)=0.d0
      MN2(7,1)=0.d0
      MN2(6,2)=0.d0
      MN2(1,3)=dsqrt(g1q/2d0)*(M1r*vuq*DDCOS(PhiM1+phiH/2d0)
     .                           +mur*vdq*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(3,1)=dsqrt(g1q/2d0)*(M1r*vuq*DDCOS(PhiM1)
     .                           +mur*vdq*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(6,8)=dsqrt(g1q/2d0)*(M1r*vuq*DDCOS(PhiM1)
     .                           +mur*vdq*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(8,6)=dsqrt(g1q/2d0)*(M1r*vuq*DDCOS(PhiM1)
     .                           +mur*vdq*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(1,8)=dsqrt(g1q/2d0)*(-M1r*vuq*DDSIN(PhiM1+phiH/2d0)
     .                           +mur*vdq*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(3,6)=dsqrt(g1q/2d0)*(M1r*vuq*DDSIN(PhiM1+phiH/2d0)
     .                           -mur*vdq*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(8,1)=dsqrt(g1q/2d0)*(-M1r*vuq*DDSIN(PhiM1+phiH/2d0)
     .                           +mur*vdq*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(6,3)=dsqrt(g1q/2d0)*(M1r*vuq*DDSIN(PhiM1+phiH/2d0)
     .                           -mur*vdq*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(1,4)=-dsqrt(g1q/2d0)*(M1r*vdq*DDCOS(PhiM1+phiH/2d0)
     .                           +mur*vuq*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(4,1)=-dsqrt(g1q/2d0)*(M1r*vdq*DDCOS(PhiM1+phiH/2d0)
     .                           +mur*vuq*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(6,9)=-dsqrt(g1q/2d0)*(M1r*vdq*DDCOS(PhiM1+phiH/2d0)
     .                           +mur*vuq*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(9,6)=-dsqrt(g1q/2d0)*(M1r*vdq*DDCOS(PhiM1+phiH/2d0)
     .                           +mur*vuq*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(1,9)=dsqrt(g1q/2d0)*(M1r*vdq*DDSIN(PhiM1+phiH/2d0)
     .                           -mur*vuq*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(9,1)=dsqrt(g1q/2d0)*(M1r*vdq*DDSIN(PhiM1+phiH/2d0)
     .                           -mur*vuq*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(4,6)=dsqrt(g1q/2d0)*(-M1r*vdq*DDSIN(PhiM1+phiH/2d0)
     .                           +mur*vuq*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(6,4)=dsqrt(g1q/2d0)*(-M1r*vdq*DDSIN(PhiM1+phiH/2d0)
     .                           +mur*vuq*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(5,1)=0.d0
      MN2(1,5)=0.d0
      MN2(6,10)=0.d0
      MN2(10,6)=0.d0
      MN2(10,1)=0.d0
      MN2(1,10)=0.d0
      MN2(10,1)=0.d0
      MN2(1,10)=0.d0
      MN2(6,5)=0.d0
      MN2(5,6)=0.d0
      MN2(2,2)=M2r**2+g2q/2d0*(vuq**2+vdq**2)
      MN2(7,7)=M2r**2+g2q/2d0*(vuq**2+vdq**2)
      MN2(2,7)=0.d0
      MN2(7,2)=0.d0
      MN2(2,3)=-dsqrt(g2q/2d0)*(M2r*vuq*DDCOS(PhiM2+phiH/2d0)
     .                           +mur*vd*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(3,2)=-dsqrt(g2q/2d0)*(M2r*vuq*DDCOS(PhiM2+phiH/2d0)
     .                           +mur*vd*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(7,8)=-dsqrt(g2q/2d0)*(M2r*vuq*DDCOS(PhiM2+phiH/2d0)
     .                           +mur*vd*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(8,7)=-dsqrt(g2q/2d0)*(M2r*vuq*DDCOS(PhiM2+phiH/2d0)
     .                           +mur*vd*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(2,8)=dsqrt(g2q/2d0)*(M2r*vuq*DDSIN(PhiM2+phiH/2d0)
     .                           -mur*vd*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(8,2)=dsqrt(g2q/2d0)*(M2r*vuq*DDSIN(PhiM2+phiH/2d0)
     .                           -mur*vd*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(3,7)=dsqrt(g2q/2d0)*(-M2r*vuq*DDSIN(PhiM2+phiH/2d0)
     .                           +mur*vd*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(7,3)=dsqrt(g2q/2d0)*(-M2r*vuq*DDSIN(PhiM2+phiH/2d0)
     .                           +mur*vd*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(2,4)=dsqrt(g2q/2d0)*(M2r*vdq*DDCOS(PhiM2+phiH/2d0)
     .                           +mur*vu*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(4,2)=dsqrt(g2q/2d0)*(M2r*vdq*DDCOS(PhiM2+phiH/2d0)
     .                           +mur*vu*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(7,9)=dsqrt(g2q/2d0)*(M2r*vdq*DDCOS(PhiM2+phiH/2d0)
     .                           +mur*vu*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(9,7)=dsqrt(g2q/2d0)*(M2r*vdq*DDCOS(PhiM2+phiH/2d0)
     .                           +mur*vu*DDCOS(Phi01+phiS+phiH/2d0))
      MN2(2,9)=dsqrt(g2q/2d0)*(-M2r*vdq*DDSIN(PhiM2+phiH/2d0)
     .                           +mur*vu*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(9,2)=dsqrt(g2q/2d0)*(-M2r*vdq*DDSIN(PhiM2+phiH/2d0)
     .                           +mur*vu*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(7,4)=dsqrt(g2q/2d0)*(M2r*vdq*DDSIN(PhiM2+phiH/2d0)
     .                           -mur*vu*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(4,7)=dsqrt(g2q/2d0)*(M2r*vdq*DDSIN(PhiM2+phiH/2d0)
     .                           -mur*vu*DDSIN(Phi01+phiS+phiH/2d0))
      MN2(2,5)=0.d0
      MN2(5,2)=0.d0
      MN2(7,10)=0.d0
      MN2(10,7)=0.d0
      MN2(2,10)=0.d0
      MN2(10,2)=0.d0
      MN2(5,7)=0.d0
      MN2(7,5)=0.d0
      MN2(3,3)=mur**2+l**2*vdq**2+(g1q+g2q)/2d0*vuq**2
      MN2(8,8)=mur**2+l**2*vdq**2+(g1q+g2q)/2d0*vuq**2
      MN2(3,8)=0.d0
      MN2(8,3)=0.d0
      MN2(3,4)=(l**2-(g1q+g2q)/2d0)*vuq*vdq
      MN2(4,3)=(l**2-(g1q+g2q)/2d0)*vuq*vdq
      MN2(8,9)=(l**2-(g1q+g2q)/2d0)*vuq*vdq
      MN2(9,8)=(l**2-(g1q+g2q)/2d0)*vuq*vdq
      MN2(3,9)=0.d0
      MN2(9,3)=0.d0
      MN2(4,8)=0.d0
      MN2(8,4)=0.d0
      MN2(3,5)=l*mur*vuq*DDCOS(phiS-phiH/2d0)
     .        -l*vdq*(mupsi*DDCOS(Phi01-phip+phiH/2d0)
     .               +ks2si*DDCOS(Phi01-Phi02+phiH/2d0-phiS))
      MN2(5,3)=l*mur*vuq*DDCOS(phiS-phiH/2d0)
     .        -l*vdq*(mupsi*DDCOS(Phi01-phip+phiH/2d0)
     .               +ks2si*DDCOS(Phi01-Phi02+phiH/2d0-phiS))
      MN2(8,10)=l*mur*vuq*DDCOS(phiS-phiH/2d0)
     .        -l*vdq*(mupsi*DDCOS(Phi01-phip+phiH/2d0)
     .               +ks2si*DDCOS(Phi01-Phi02+phiH/2d0-phiS))
      MN2(10,8)=l*mur*vuq*DDCOS(phiS-phiH/2d0)
     .        -l*vdq*(mupsi*DDCOS(Phi01-phip+phiH/2d0)
     .               +ks2si*DDCOS(Phi01-Phi02+phiH/2d0-phiS))
      MN2(3,10)=-l*mur*vuq*DDSIN(phiS-phiH/2d0)
     .          +l*vdq*(mupsi*DDSIN(Phi01-phip+phiH/2d0)
     .                 +ks2si*DDSIN(Phi01-Phi02+phiH/2d0-phiS))
      MN2(10,3)=-l*mur*vuq*DDSIN(phiS-phiH/2d0)
     .          +l*vdq*(mupsi*DDSIN(Phi01-phip+phiH/2d0)
     .                 +ks2si*DDSIN(Phi01-Phi02+phiH/2d0-phiS))
      MN2(5,8)=l*mur*vuq*DDSIN(phiS-phiH/2d0)
     .         -l*vdq*(mupsi*DDSIN(Phi01-phip+phiH/2d0)
     .                 +ks2si*DDSIN(Phi01-Phi02+phiH/2d0-phiS))
      MN2(8,5)=l*mur*vuq*DDSIN(phiS-phiH/2d0)
     .         -l*vdq*(mupsi*DDSIN(Phi01-phip+phiH/2d0)
     .                 +ks2si*DDSIN(Phi01-Phi02+phiH/2d0-phiS))
      MN2(4,4)=mur**2+l**2*vuq**2+(g1q+g2q)/2d0*vdq**2
      MN2(9,9)=mur**2+l**2*vuq**2+(g1q+g2q)/2d0*vdq**2
      MN2(4,9)=0.d0
      MN2(9,4)=0.d0
      MN2(4,5)=l*mur*vdq*DDCOS(phiS-phiH/2d0)
     .        -l*vuq*(mupsi*DDCOS(Phi01-phip+phiH/2d0)
     .               +ks2si*DDCOS(Phi01-Phi02+phiH/2d0-phiS))
      MN2(5,4)=l*mur*vdq*DDCOS(phiS-phiH/2d0)
     .        -l*vuq*(mupsi*DDCOS(Phi01-phip+phiH/2d0)
     .               +ks2si*DDCOS(Phi01-Phi02+phiH/2d0-phiS))
      MN2(9,10)=l*mur*vdq*DDCOS(phiS-phiH/2d0)
     .        -l*vuq*(mupsi*DDCOS(Phi01-phip+phiH/2d0)
     .               +ks2si*DDCOS(Phi01-Phi02+phiH/2d0-phiS))
      MN2(10,9)=l*mur*vdq*DDCOS(phiS-phiH/2d0)
     .        -l*vuq*(mupsi*DDCOS(Phi01-phip+phiH/2d0)
     .               +ks2si*DDCOS(Phi01-Phi02+phiH/2d0-phiS))
      MN2(4,10)=-l*mur*vdq*DDSIN(phiS-phiH/2d0)
     .          +l*vuq*(mupsi*DDSIN(Phi01-phip+phiH/2d0)
     .               +ks2si*DDSIN(Phi01-Phi02+phiH/2d0-phiS))
      MN2(10,4)=-l*mur*vdq*DDSIN(phiS-phiH/2d0)
     .          +l*vuq*(mupsi*DDSIN(Phi01-phip+phiH/2d0)
     .               +ks2si*DDSIN(Phi01-Phi02+phiH/2d0-phiS))
      MN2(5,9)=l*mur*vdq*DDSIN(phiS-phiH/2d0)
     .         -l*vuq*(mupsi*DDSIN(Phi01-phip+phiH/2d0)
     .               +ks2si*DDSIN(Phi01-Phi02+phiH/2d0-phiS))
      MN2(9,5)=l*mur*vdq*DDSIN(phiS-phiH/2d0)
     .         -l*vuq*(mupsi*DDSIN(Phi01-phip+phiH/2d0)
     .               +ks2si*DDSIN(Phi01-Phi02+phiH/2d0-phiS))
      MN2(5,5)=mupsi**2+ks2si**2+2d0*mupsi*ks2si*DDCOS(phip-phi02-phiS)
     .        +l**2*(vuq**2+vdq**2)
      MN2(10,10)=mupsi**2+ks2si**2+2d0*mupsi*ks2si
     .        *DDCOS(phip-phi02-phiS)+l**2*(vuq**2+vdq**2)
      MN2(5,10)=0.d0
      MN2(10,5)=0.d0
      CALL DIAGN(10,MN2,VALPN,VECPN,1d-10)
      CALL SORTNA(10,VALPN,VECPN)

      Taux=MHuS+MHdS+2d0*muq**2+g2q/2d0*(vuq**2+vdq**2)
      Daux=MHuS-MHdS+g1q/2d0*(vuq**2-vdq**2)
      aux=M3HQ**2+(l*XIFQ)**2+2d0*M3HQ*l*XIFQ*DDCOS(phi3q+phiF-phi01)
     .   +(MUPQ*muq)**2+2d0*MUPQ*muq*(M3HQ*DDCOS(phi3q+phiP-phi01+phiS)
     .                                    +l*XIFQ*DDCOS(phiP-phiF+phiS))
     .   +muq**2*(Alcos1**2+IAl**2+(k/l*muq)**2+2d0*k/l*muq
     .       *(Alcos1*DDCOS(phi0-3d0*phiS)+IAl*DDSIN(phi0-3d0*phiS)))
     .   +(l**2-g2q/2d0)**2*(vuq*vdq)**2
     . -2d0*(l**2-g2q/2d0)*muq*vuq*vdq*(Alcos1*DDCOS(phiH+phiS)
     .       -IAl*DDSIN(phiH+phiS)+k/l*muq*DDCOS(phi0+phiH-2d0*phiS))
     . +2d0*M3HQ*(muq*(Alcos1*DDCOS(phi3q-phiS)+IAL*DDSIN(phi3q-phiS))
     .           +k/l*muq**2*DDCOS(phi3q-phi0+2d0*phiS)
     .           -(l**2-g2q/2d0)*vuq*vdq*DDCOS(phi3q+phiH))
     . +2d0*l*XIFQ*(muq*(Alcos1*DDCOS(phi01-phiF-phiS)
     .                       +IAL*DDSIN(phi01-phiF-phiS))
     .           +k/l*muq**2*DDCOS(phi02-phiF+2d0*phiS)
     .           -(l**2-g2q/2d0)*vuq*vdq*DDCOS(phi01-phiF+phiH))
     . +2d0*MUPQ*muq*(muq*(Alcos1*DDCOS(phi01-phiP-2d0*phiS)
     .                       +IAL*DDSIN(phi01-phiP-2d0*phiS))
     .           +k/l*muq**2*DDCOS(phi02-phiP+3d0*phiS)
     .           -(l**2-g2q/2d0)*vuq*vdq*DDCOS(phi01-phiP+phiS+phiH))
      MHD1aux=(Taux-dsqrt(Daux**2+4d0*aux))/2d0
      MHD2aux=(Taux+dsqrt(Daux**2+4d0*aux))/2d0

      MH2(1,1)=MHuS+muq**2+(l**2-(g1q+g2q)/4d0)*vdq**2
     .        +(g1q+g2q)/4d0*vuq**2*(1d0+DDCOS(phiH/2d0)**2)
      MH2(1,2)=-muq*(Alcos1*DDCOS(phiS)-IAl*DDSIN(phiS))
     .        -k/l*muq**2*DDCOS(phi0-2d0*phiS)
     .        +2d0*(l**2-(g1q+g2q)/4d0)*vuq*vdq*DDCOS(phiH/2d0)**2
     .        -M3HQ*DDCOS(phi3q)-l*XIFQ*DDCOS(phi01-phiF)
     .        -muq*MUPQ*DDCOS(phi01-phiP-phiS)
      MH2(2,1)=MH2(1,2)
      MH2(1,3)=-l*vdq*(Alcos1*DDCOS(phiH/2d0)-IAl*DDSIN(phiH/2d0)
     .                +2d0*k/l*muq*DDCOS(phi0-phiS+phiH/2d0)
     .                +MUPQ*DDCOS(phi01-phiP+phiH/2))
     .         +2d0*l*muq*vuq*DDCOS(phiS)*DDSIN(phiH/2d0)
      MH2(3,1)=MH2(1,3)
      MH2(2,2)=MHdS+muq**2+(l**2-(g1q+g2q)/4d0)*vuq**2
     .        +(g1q+g2q)/4d0*vdq**2*(1d0+DDCOS(phiH/2d0)**2)
      MH2(2,3)=-l*vuq*(Alcos1*DDCOS(phiH/2d0)-IAl*DDSIN(phiH/2d0)
     .                +2d0*k/l*muq*DDCOS(phi0-phiS+phiH/2d0)
     .                +MUPQ*DDCOS(phi01-phiP+phiH/2d0))
     .         +2d0*l*muq*vdq*DDCOS(phiS)*DDSIN(phiH/2d0)
      MH2(3,2)=MH2(2,3)
      MH2(3,3)=MSS+l**2*(vuq**2+vdq**2)+MUPQ**2
     .        +2d0*k/l*muq*(Akcos2*DDCOS(phiS)-IAk*DDSIN(phiS))
     .        +2d0*(k/l*muq)**2*(1d0+1d0*DDCOS(phiS)**2)
     .        -2d0*l*k*vuq*vdq*DDCOS(phi0+phiH)
     .        +MSPQ*DDCOS(phiSPq)+2d0*k*XIFQ*DDCOS(phi02-phiF)
     .        +2d0*k/l*muq*MUPQ*(DDCOS(phi02-phiP+phiS)
     .                          +2d0*DDCOS(phi02-phiP)*DDCOS(phiS))
      MH2(1,4)=(g1q+g2q)/4d0*vuq**2*DDSIN(phiH)
      MH2(4,1)=MH2(1,4)
      MH2(2,4)=muq*(IAl*DDCOS(phiS)+Alcos1*DDSIN(phiS)
     .              +k/l*muq*DDSIN(phi0-2d0*phiS))
     .        +(l**2-(g1q+g2q)/4d0)*vuq*vdq*DDSIN(phiH)
     .        +M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
     .        +muq*MUPQ*DDSIN(phi01-phiP-phiS)
      MH2(4,2)=MH2(2,4)
      MH2(3,4)=l*vdq*(IAl*DDCOS(phiH/2d0)+Alcos1*DDSIN(phiH/2d0)
     .               +2d0*k/l*muq*DDSIN(phi0-phiS+phiH/2d0)
     .               +MUPQ*DDSIN(phi01-phiP+phiH/2d0))
     .        +2d0*l*muq*DDCOS(phiS)*DDSIN(phiH/2d0)
      MH2(4,3)=MH2(3,4)
      MH2(4,4)=MHuS+muq**2+(l**2-(g1q+g2q)/4d0)*vdq**2
     .        +(g1q+g2q)/4d0*vuq**2*(1d0+2d0*DDSIN(phiH/2d0)**2)
      MH2(1,5)=muq*(IAl*DDCOS(phiS)+Alcos1*DDSIN(phiS)
     .             +k/l*muq*DDSIN(phi0-2d0*phiS))
     .       +(l**2-(g1q+g2q)/4d0)*vuq*vdq*DDSIN(phiH)
     .       +M3HQ*DDSIN(phi3q)+l*XIFQ*DDSIN(phi01-phiF)
     .       +muq*MUPQ*DDSIN(phi01-phiP-phiS)
      MH2(5,1)=MH2(1,5)
      MH2(2,5)=(g1q+g2q)/4d0*vdq**2*DDSIN(phiH)
      MH2(5,2)=MH2(2,5)
      MH2(3,5)=l*vuq*(IAl*DDCOS(phiH/2d0)+Alcos1*DDSIN(phiH/2d0)
     .               +2d0*k/l*muq*DDSIN(phi0-phiS+phiH/2d0)
     .               +MUPQ*DDSIN(phi01-phiP+phiH/2d0))
     .         +2d0*l*muq*vdq*DDCOS(phiS)*DDSIN(phiH/2d0)
      MH2(5,3)=MH2(3,5)
      MH2(4,5)=muq*(Alcos1*DDCOS(phiS)-IAl*DDSIN(phiS)
     .             +k/l*muq*DDCOS(phi0-2d0*phiS))
     .        +2d0*(l**2-(g1q+g2q)/4d0)*vuq*vdq*DDSIN(phiH/2d0)**2
     .        +M3HQ*DDCOS(phi3q)+l*XIFQ*DDCOS(phi01-phiF)
     .        +muq*MUPQ*DDCOS(phi01-phiP-phiS)
      MH2(5,4)=MH2(4,5)
      MH2(5,5)=MHdS+muq**2+(l**2-(g1q+g2q)/4d0)*vuq**2
     .        +(g1q+g2q)/4d0*vdq**2*(1d0+2d0*DDSIN(phiH/2d0)**2)
      MH2(1,6)=l*vdq*(IAl*DDCOS(phiH/2d0)+Alcos1*DDSIN(phiH/2d0)
     .                -2d0*k/l*muq*DDSIN(phi0-phiS+phiH/2d0)
     .                +MUPQ*DDSIN(phi01-phiP+phiH/2d0))
     .         +2d0*l*muq*vuq*DDSIN(phiS)*DDCOS(phiH/2d0)
      MH2(6,1)=MH2(1,6)
      MH2(2,6)=l*vuq*(IAl*DDCOS(phiH/2d0)+Alcos1*DDSIN(phiH/2d0)
     .               -2d0*k/l*muq*DDSIN(phi0-phiS+phiH/2d0)
     .               +MUPQ*DDSIN(phi01-phiP+phiH/2d0))
     .         +2d0*l*muq*vdq*DDSIN(phiS)*DDCOS(phiH/2d0)
      MH2(6,2)=MH2(2,6)
      MH2(3,6)=-2d0*k/l*muq*(IAk*DDCOS(phiS)+Akcos2*DDSIN(phiS)
     .                      -k/l*muq*DDSIN(2d0*phiS))
     .   -2d0*k*l*vuq*vdq*DDSIN(phi0+phiH)
     .   -MSPQ*DDSIN(phiSPq)-2d0*k*XIFQ*DDSIN(phi02-phiF)
     .   +2d0*k/l*muq*MUPQ*(-DDSIN(phi02-phiP)*DDCOS(phiS)
     .                        +DDCOS(phi02-phiP)*DDSIN(phiS))
      MH2(6,3)=MH2(3,6)
      MH2(4,6)=l*vdq*(Alcos1*DDCOS(phiH/2d0)-IAl*DDSIN(phiH/2d0)
     .               -2d0*k/l*muq*DDCOS(phi0-phiS+phiH/2d0)
     .               +MUPQ*DDCOS(phi01-phiP+phiH/2d0))
     .        +2d0*l*muq*vuq*DDSIN(phiS)*DDCOS(phiH/2d0)
      MH2(6,4)=MH2(4,6)
      MH2(5,6)=l*vuq*(Alcos1*DDCOS(phiH/2d0)-IAl*DDSIN(phiH/2d0)
     .               -2d0*k/l*muq*DDCOS(phi0-phiS+phiH/2d0)
     .               +MUPQ*DDCOS(phi01-phiP+phiH/2d0))
     .        +2d0*l*muq*vdq*DDSIN(phiS)*DDCOS(phiH/2d0)
      MH2(6,5)=MH2(4,5)
      MH2(6,6)=MSS-2d0*k/l*muq*(Akcos2*DDCOS(phiS)-IAk*DDSIN(phiS)
     .                         +k/l*muq*(1d0+2d0*DDSIN(phiS)**2))
     .      +l**2*(vuq**2+vdq**2)+2d0*l*k*vuq*vdq*DDCOS(phi0+phiH)
     .      +MUPQ**2-MSPQ*DDCOS(phiSPq)-2d0*k*XIFQ*DDCOS(phi02-phiF)
     .      +2d0*k/l*muq*MUPQ*(DDCOS(phi02-phiP+phiS)
     .                        -2d0*DDSIN(phi02-phiP)*DDSIN(phiS))
      CALL DIAGN(6,MH2,VALPH,VECPH,1d-10)

      Vphi=min(Vphi,MHuS*vuq**2+MHdS*vdq**2+MSS*(muq/l)**2
     . -2d0*muq*vuq*vdq*(Alcos1*DDCOS(phiS+phiH)-IAl*DDSIN(phiS+phiH))
     . +2d0/3d0*k*(muq/l)**3*(Akcos2*DDCOS(3*PhiS)-IAk*DDSIN(3d0*phiS))
     . +muq**2*(vuq**2+vdq**2)+(l*vuq*vdq)**2
     . -2d0*k/l*muq**2*vuq*vdq*DDCOS(phi0+phiH-2d0*phiS)
     . +k**2*(muq/l)**4+(G1Q+G2Q)/8d0*(vuq**2-vdq**2)**2
     . -2d0*(M3HQ*DDCOS(phi3q+phiH)+l*XIFQ*DDCOS(phi01-phiF+phiH)
     .                 +MUPQ*muq*DDCOS(phi01-phiP+phiH-phiS))*vuq*vdq
     . +2d0*(XISQ*DDCOS(phiSq)*DDCOS(phiS)-IXIS*DDSIN(PhiS)
     .                   +XIFQ*MUPQ*DDCOS(phiP-phiF+phiS))*muq/l
     . +2d0*(MSPQ/2d0*DDCOS(phiSPq+2d0*phiS)
     .                    +k*XIFQ*DDCOS(phi02-phiF+2d0*phiS))*(muq/l)**2
     . +2d0*k*MUPQ*DDCOS(phi02-phiP+phiS)*(muq/l)**3
     . +(MUPQ*muq/l)**2
     . -(3d0*mtopq**4*(dlog(mtopq**2/QSTSB)-1.5d0)
     .  +3d0*mbotq**4*(dlog(mbotq**2/QSTSB)-1.5d0)
     .  +(mtau*vdq/vd)**4*(dlog(mtau**2/QSTSB)-1.5d0))/16d0/Pi**2
     . +3d0*(2d0*MW2**2*(dlog(MW2/QSTSB)-1.5d0)
     .             +MZ2*(dlog(MZ2/QSTSB)-1.5d0))/64d0/Pi**2
     . +(3d0*MSTaux1**2*(dlog(Max(MSTaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSTaux2**2*(dlog(Max(MSTaux2,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux1**2*(dlog(Max(MSBaux1,MZ2)/QSTSB)-1.5d0)
     .     +3d0*MSBaux2**2*(dlog(Max(MSBaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux1**2*(dlog(Max(MSLaux1,MZ2)/QSTSB)-1.5d0)
     .     +MSLaux2**2*(dlog(Max(MSLaux2,MZ2)/QSTSB)-1.5d0)
     .     +MSNTaux**2*(dlog(Max(MSNTaux,MZ2)/QSTSB)-1.5d0)
     .     +6d0*MSU2(1)**2*(dlog(MSU2(1)/QSTSB)-1.5d0)
     .     +6d0*MSU2(2)**2*(dlog(MSU2(2)/QSTSB)-1.5d0)
     .     +6d0*MSD2(1)**2*(dlog(MSD2(1)/QSTSB)-1.5d0)
     .     +6d0*MSD2(2)**2*(dlog(MSD2(2)/QSTSB)-1.5d0)
     .     +2d0*MSE2(1)**2*(dlog(MSE2(1)/QSTSB)-1.5d0)
     .     +2d0*MSE2(2)**2*(dlog(MSE2(2)/QSTSB)-1.5d0)
     .     +2d0*MSNE2**2*(dlog(MSNE2/QSTSB)-1.5d0))/32d0/Pi**2
     .    +3d0/512d0/Pi**4*Ytq**4*((dlog(QSTSB/mtopq**2))**2
     .      *(64d0*Pi*ALSQ+4d0/3d0*g1q-3d0*sinb**2*Ytq**2
     .       +3d0*cosb**2*Ybq**2)+
     .       ((dlog(MA2/mtopq**2))**2-(dlog(QSTSB/mtopq**2))**2)
     .       *(3d0*cosb**2*Ytq**2+(3d0*cosb**2+1d0)*Ybq**2))*vuq**4
     .    +3d0/512d0/Pi**4*Ybq**4*(dlog(QSTSB/mtopq**2)**2
     .      *(64d0*Pi*ALSQ-2d0/3d0*g1q+3d0*sinb**2*Ytq**2
     .       -3d0*cosb**2*Ybq**2)+
     .       (dlog(MA2/mtopq**2)**2-dlog(QSTSB/mtopq**2)**2)
     .       *(3d0*sinb**2*Ybq**2+(3d0*sinb**2+1d0)*Ytq**2))*vdq**4
     .    -(MCH1aux**2*(dlog(Max(MCH1aux,MZ2)/QSTSB)-1.5d0)
     .     +MCH2aux**2*(dlog(Max(MCH2aux,MZ2)/QSTSB)-1.5d0))/16d0/Pi**2
     .    -(VALPN(1)**2*(dlog(Max(VALPN(1),1d0)/QSTSB)-1.5d0)
     .     +VALPN(3)**2*(dlog(Max(VALPN(3),1d0)/QSTSB)-1.5d0)
     .     +VALPN(5)**2*(dlog(Max(VALPN(5),1d0)/QSTSB)-1.5d0)
     .     +VALPN(7)**2*(dlog(Max(VALPN(7),1d0)/QSTSB)-1.5d0)
     .   +VALPN(9)**2*(dlog(Max(VALPN(9),1d0)/QSTSB)-1.5d0))/32d0/Pi**2
     .     +(MHD1aux**2*(dlog(Max(MHD1aux,MW2)/QSTSB)-1.5d0)
     .     +MHD2aux**2*(dlog(Max(MHD2aux,MW2)/QSTSB)-1.5d0))/32d0/Pi**2
     .    +(VALPH(1)**2*(dlog(Max(MZ2,VALPH(1))/QSTSB)-1.5d0)
     .     +VALPH(2)**2*(dlog(Max(MZ2,VALPH(2))/QSTSB)-1.5d0)
     .     +VALPH(3)**2*(dlog(Max(MZ2,VALPH(3))/QSTSB)-1.5d0)
     .     +VALPH(4)**2*(dlog(Max(MZ2,VALPH(4))/QSTSB)-1.5d0)
     .     +VALPH(5)**2*(dlog(Max(MZ2,VALPH(5))/QSTSB)-1.5d0)
     .   +VALPH(6)**2*(dlog(Max(MZ2,VALPH(6))/QSTSB)-1.5d0))/64d0/Pi**2)

      ENDDO
      ENDDO

      IF(Vphi.lt.VEW)PROB(28)=Max(PROB(28),
     .                     DDIM(-1d-2,(Vphi-VEW)/Max(dabs(VEW),1d-10)))

c      Naturalness of the doublet squared masses

      PROB(29)=DDIM(MAX(DABS(MHuS),DABS(MHdS))/Q2,10d0)

 630  RETURN
      END


text:
subroutine updveg3d(mmax      ,nmax      ,kmax      ,sig       ,thick     , &
                   & dps       ,kfs       ,s1        ,u1        ,v1        , &
                   & diapl     ,rnpl      ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: updveg3d.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/timedep/updveg3d.f90 $
!!--description-----------------------------------------------------------------
! NONE
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use mathconsts
    use globaldata
    use dfparall

    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    real(fp)                             , pointer :: amiss
    integer, dimension(:,:)              , pointer :: planttype
    real(fp), dimension(:,:)             , pointer :: nplants
    type (dpm_vegetation), dimension(:)  , pointer :: vegs
    integer                              , pointer :: lundia
    real(fp)                             , pointer :: ag
    real(fp)                             , pointer :: rhow
!
! Constants
!
    integer , parameter :: numiter = 30        ! number of iterations for calculating phi
    real(fp), parameter :: px      = 0.15_fp
    real(fp), parameter :: pm      = 0.001_fp
!
! Global variables
!
    integer                                                            , intent(in)  :: kmax
    integer                                                            , intent(in)  :: mmax
    integer                                                            , intent(in)  :: nmax
    integer, dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)       , intent(in)  :: kfs
    real(fp), dimension(kmax)                                          , intent(in)  :: sig
    real(fp), dimension(kmax)                                          , intent(in)  :: thick
    real(prec), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)    , intent(in)  :: dps
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)      , intent(in)  :: s1
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax), intent(in)  :: u1
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax), intent(in)  :: v1
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax), intent(out) :: diapl
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax), intent(out) :: rnpl
!
! Local variables
!
    integer           :: i
    integer           :: k
    integer           :: l
    integer           :: m
    integer           :: n
    integer           :: vegtype
    real(fp)          :: a
    real(fp)          :: apln
    real(fp)          :: aplo
    real(fp)          :: aplt
    real(fp)          :: arm
    real(fp)          :: aspn
    real(fp)          :: aspo
    real(fp)          :: aspt
    real(fp)          :: b
    real(fp)          :: bendm
    real(fp)          :: buoym
    real(fp)          :: cdn
    real(fp)          :: cos_angle
    real(fp)          :: dian
    real(fp)          :: diao
    real(fp)          :: dz
    real(fp)          :: dzf
    real(fp)          :: el
    real(fp)          :: fb
    real(fp)          :: fbuoy
    real(fp)          :: frac
    real(fp)          :: height
    real(fp)          :: pdens
    real(fp)          :: ph
    real(fp)          :: phi         ! angle between vegatation and z-axis
    real(fp)          :: phi_old
    real(fp)          :: phl
    real(fp)          :: pi4         ! pi/4
    real(fp)          :: rhop
    real(fp)          :: rlbot
    real(fp)          :: rltop
    real(fp)          :: rnn
    real(fp)          :: rno
    real(fp)          :: sin_angle
    real(fp)          :: tt
    real(fp)          :: u
    real(fp)          :: v
    real(fp)          :: vegcount
    real(fp)          :: vmag2
    real(fp)          :: zfl
    real(fp)          :: zflb
    real(fp)          :: zflow
    real(fp)          :: zflt
    real(fp)          :: zplb
    real(fp)          :: zplt
    character(80)     :: message
!
!! executable statements -------------------------------------------------------
!
    amiss      => gdp%gdconst%amiss
    planttype  => gdp%gdveg3d%planttype
    nplants    => gdp%gdveg3d%nplants
    vegs       => gdp%gdveg3d%vegs
    lundia     => gdp%gdinout%lundia
    ag         => gdp%gdphysco%ag
    rhow       => gdp%gdphysco%rhow
    pi4 = pi * 0.25_fp
    !
    ! Initialization
    !
    do k = 1, kmax
       do n = 1, nmax
          do m = 1, mmax
             diapl(n, m, k) = 0.0_fp
             rnpl(n, m, k)  = 0.0_fp
          enddo
       enddo
    enddo
    do m = 1, mmax
       do n = 1, nmax
          if (abs(kfs(n, m)) == 1) then
             if (planttype(n,m) /= 0) then
                vegtype  = planttype(n,m)
                vegcount = nplants(n,m)
                height   = s1(n, m) + real(dps(n, m),fp)
                !
                ! Calculate angle of plant
                ! Assumptions:
                ! - The vegetation is fully stiff
                ! - The vegetation is "tied" to the bottom
                !
                if (comparereal(vegs(vegtype)%rho(1),amiss) == 0) then
                   phi = 0.0_fp
                else
                   phi = 0.3_fp
                   do i = 1,numiter
                      phi_old   = phi
                      cos_angle = cos(phi)
                      sin_angle = sin(phi)
                      bendm     = 0.0_fp
                      buoym     = 0.0_fp
                      do k  = 1,kmax
                         u     = (u1(n,m,k) + u1(n,m-1,k)) / 2.0_fp
                         v     = (v1(n,m,k) + v1(n-1,m,k)) / 2.0_fp
                         vmag2 = u*u+v*v
                         zflow  = (1.0_fp+sig(k)) * height       !center layer
                         dzf    = thick(k)   * height            !thickness layer
                         zflt   = zflow + dzf/2.0_fp             !upper boundary layer
                         zflb   = zflow - dzf/2.0_fp + 0.0001_fp !lower boundary layer
                         do l = 1,vegs(vegtype)%nvps-1
                            zplt = vegs(vegtype)%z(l+1)* cos_angle       !upper side plantelement
                            zplb = vegs(vegtype)%z(l)  * cos_angle       !lower side plantelement
                            if (      comparereal(zflb,zplb) /= -1 &
                              & .and. comparereal(zflb,zplt) /=  1) then
                               a     =   (zflow                - vegs(vegtype)%z(l)) &
                                     & / (vegs(vegtype)%z(l+1) - vegs(vegtype)%z(l))
                               b     =  1.0_fp - a
                               cdn   = (b*vegs(vegtype)%cdcoef(l) + a*vegs(vegtype)%cdcoef(l+1))
                               dian  = (b*vegs(vegtype)%dia   (l) + a*vegs(vegtype)%dia   (l+1))*cdn
                               rnn   = (b*vegs(vegtype)%nstem (l) + a*vegs(vegtype)%nstem (l+1))*vegcount
                               fb    = 0.5_fp * cdn * dian * rnn * vmag2 * rhow
                               frac  = (zplt - max(zplb,zflb)) / dzf  !fraction plant in layer
                               frac  = min(1.0_fp,frac)
                               fb    = fb * thick(k) * height * frac
                               arm   = zflow - 0.5_fp*thick(k)*height*(1.0_fp-frac)
                               bendm = bendm + fb*arm
                            endif
                         enddo
                      enddo
                      arm   = 0.0_fp
                      buoym = 0.0_fp
                      do l  = 1,vegs(vegtype)%nvps-1  ! then upwelling moment buoym
                         dian  = 0.5_fp * (vegs(vegtype)%dia  (l) + vegs(vegtype)%dia  (l+1))
                         rnn   = 0.5_fp * (vegs(vegtype)%nstem(l) + vegs(vegtype)%nstem(l+1))
                         rhop  = 0.5_fp * (vegs(vegtype)%rho  (l) + vegs(vegtype)%rho  (l+1))
                         apln  = dian * dian * rnn * pi4
                         rltop = min(vegs(vegtype)%z(l+1) , height/cos_angle)
                         rlbot = min(vegs(vegtype)%z(l  ) , height/cos_angle)
                         el    = rltop  - rlbot
                         fbuoy = (rhow  - rhop) * ag * apln * el
                         arm   = rltop  - 0.5_fp*el
                         buoym = buoym + fbuoy*arm*sin_angle
                      enddo
                      if (comparereal(bendm,0.0_fp) /= 0) then
                         ph  = (bendm - buoym) / (bendm + buoym)
                         tt  = real(i-1,fp) / real(numiter,fp)
                         phl = tt*pm + (1-tt)*px
                         ph  = max( -phl , min(phl,ph) )
                      else
                         ph = 0.0_fp
                      endif
                      phi = phi + ph
                      phi = min( pi*0.5_fp , max(0.0_fp,phi) )
                      !
                      ! Jump out of iteration loop when phi has not changed much
                      !
                      if (comparereal(abs(phi-phi_old),0.01_fp) == -1) exit
                   enddo
                endif
                cos_angle = cos(phi)
                do k = 1, kmax
                   zflow = (1.0_fp+sig(k)) * height
                   dzf   = thick(k)        * height
                   zflt  = zflow + dzf/2.0_fp
                   zflb  = zflow - dzf/2.0_fp + 0.0001_fp
                   do l = 1, vegs(vegtype)%nvps - 1
                      zplt = vegs(vegtype)%z(l+1) * cos_angle
                      zplb = vegs(vegtype)%z(l)   * cos_angle
                      if (        comparereal(zflb,zplb) >  -1 &
                          & .and. comparereal(zflb,zplt) == -1  ) then
                         dz = zplt - zplb
                         if (comparereal(dz,0.0_fp) == 0) then
                            write(message,'(a,i0)') 'Identical z coordinates in specification of vegetation type ',vegtype
                            call prterr(lundia, 'U021', message)
                            call d3stop(1, gdp)
                         endif
                         zfl  = 0.5_fp * (min(zplt,zflt) + zflb)
                         a    = (zfl - zplb)/dz
                         b    = 1.0_fp - a
                         diao = diapl(n,m,k)
                         rno  = rnpl (n,m,k)
                         aspo = diao * rno                 ! aspect  old
                         aplo = diao * diao * rno * pi4    ! a plane old
                         !
                         ! Version Arjen has cdn = 1
                         !
                         cdn  = (b*vegs(vegtype)%cdcoef(l) + a*vegs(vegtype)%cdcoef(l+1))
                         !
                         ! cd is stored in dia
                         !
                         dian = (b*vegs(vegtype)%dia   (l) + a*vegs(vegtype)%dia   (l+1)) * cdn
                         rnn  = (b*vegs(vegtype)%nstem (l) + a*vegs(vegtype)%nstem (l+1)) * vegcount
                         if (       comparereal(aspo,0.0_fp) == 0 &
                             & .or. comparereal(aplo,0.0_fp) == 0  ) then
                            diapl(n,m,k) = dian
                            rnpl (n,m,k) = rnn
                         else
                            aspn         = dian * rnn               ! aspect  new
                            apln         = dian * dian * rnn * pi4  ! a plane new
                            aspt         = aspo + aspn              ! aspect  tot
                            aplt         = aplo + apln              ! a plane tot
                            diapl(n,m,k) = aplt / (aspt*pi4)
                            rnpl (n,m,k) = (aspt*aspt*pi4) / aplt
                         endif
                      endif
                   enddo
                enddo
             endif
          endif
       enddo
    enddo
    !
    ! Check plant density
    !
    do k = 1, kmax
       do n = 1, nmax
          do m = 1, mmax
             pdens = 1.0 - diapl(n, m, k)*diapl(n, m, k)*rnpl(n, m, k)*pi4
             if (comparereal(pdens,0.01_fp) == -1) then
                write(message,'(a,i0,a,i0,a)') 'plant density higher than possible for (m,n) = (', &
                    &                          m,',',n,')'
                call prterr(lundia, 'U021', message)
                call d3stop(1, gdp)
             endif
          enddo
       enddo
    enddo
end subroutine updveg3d

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      PROGRAM UAREA
C       A User's access to the ITM in the area prediction mode
C          Programming error here corrected 3/17/98
C
      CHARACTER KPSYS*6,KPENV*5,KPSTC*7,KPTAB*5
      CHARACTER LPSYS(6)*8,LPENV(5)*8,LPSTC(7)*8,LPTAB(5)*8
      DIMENSION IPSYS(6)
        EQUIVALENCE(SH1,IPSYS(1)),(SH2,IPSYS(2)),(FMHZ,IPSYS(3)),
     X    (KPOL,IPSYS(4)),(KST1,IPSYS(5)),(KST2,IPSYS(6))
      DIMENSION IPENV(5)
        EQUIVALENCE(SDH,IPENV(1)),(EN0,IPENV(2)),(KLMM,IPENV(3)),
     X    (EPS,IPENV(4)),(SGM,IPENV(5))
      DIMENSION IPSTC(7),QC(4)
        EQUIVALENCE(KVAR,IPSTC(1)),(QT,IPSTC(2)),(QL,IPSTC(3)),
     X    (QC,IPSTC(4))
      DIMENSION IPTAB(5)
        EQUIVALENCE(D0,IPTAB(1)),(D1,IPTAB(2)),(DS1,IPTAB(3)),
     X    (D2,IPTAB(4)),(DS2,IPTAB(5))
C
      COMMON/PROP/KWX,AREF,MDP,DIST,HG(2),WN,DH,ENS,GME,ZGND,
     X  HE(2),DL(2),THE(2)
         COMPLEX ZGND
      COMMON/PROPV/LVAR,SGC,MDVAR,KLIM
C
      CHARACTER LQ*1,TITLE*72,FNAME*40,FNQ*40
      DIMENSION ZC(4),XLB(4)
      DIMENSION KST(2)
C
      PARAMETER(GMA=157E-9)
      PARAMETER(DB=8.685890)
      PARAMETER(AKM=1000.)
C
      DATA LPSYS/' H1, m',' H2, m',' F, MHz','  Pol',' Sit1',' Sit2'/
        DATA SH1,SH2,FMHZ,KPOL,KST1,KST2/3.,3.,100.,86,82,82/
        DATA KPSYS/'111AAA'/
      DATA LPENV/' DH, m','  Ns',' Clim','  eps','  sgm'/
        DATA SDH,EN0,KLMM,EPS,SGM/90.,301.,67,15.,.005/
        DATA KPENV/'00A03'/
      DATA LPSTC/'  Var',' QT, %',' QL, %',
     X           '  QC1','  QC2','  QC3','  QC4'/
        DATA KVAR,QT,QL,QC/66,50.,50.,50.,90.,10.,0./
        DATA KPSTC/'A111111'/
      DATA LPTAB/' d0, km','  d1','  ds1','  d2','  ds2'/
        DATA D0,D1,DS1,D2,DS2/10.,150.,10.,500.,50./
        DATA KPTAB/'00000'/
C
      DATA TITLE/' '/,FNAME/'ITM.DAT'/,NNAME/7/,JFILE/1/
C
2000  FORMAT(A)
2002  FORMAT(2X,A,' ',$)
2102  FORMAT(8X,'At ',A6,',  F(50,',A2,') =',F6.1,' dBu')
C
      WRITE(*,2000) ' ',' '
      WRITE(*,2000) '     The Irregular Terrain Model'//
     X              ' of Radio Propagation'
      WRITE(*,2000) '     in the Area Prediction mode'
      WRITE(*,2000) ' '
      WRITE(*,2000) '     Input:'
10    CONTINUE
      WRITE(*,2000) '     System parameters--'
      WRITE(*,2000) '       Pol=H|V, Sit=R|C|V'
      CALL VCARD(6,KPSYS,LPSYS,IPSYS,*80)
      WRITE(*,2000) '     Environment parameters--'
      WRITE(*,2000) '       Clim=C|M|S|1|2...|7'
      CALL VCARD(5,KPENV,LPENV,IPENV,*80)
      WRITE(*,2000) '     Statistical parameters--'
      WRITE(*,2000) '       Var=S|A|M|B'
      CALL VCARD(7,KPSTC,LPSTC,IPSTC,*80)
      WRITE(*,2000) '     Tabulation parameters--'
      CALL VCARD(5,KPTAB,LPTAB,IPTAB,*80)
C
      WRITE(*,2002) '...needs editing (N)?'
      READ(*,2000) LQ
      IF(LQ .EQ. 'Y' .OR. LQ .EQ. 'y') GO TO 10
      WRITE(*,2002) '   A title?'
      READ(*,2000) TITLE(4:)
11    WRITE(*,2002) '   and the output file: '//FNAME(:NNAME)//'?'
      READ(*,2000) FNQ
      IF(FNQ .NE. ' ') THEN
        JFILE=1
        FNAME=' '
        NNAME=0
        DO I=1,40
          LQ=FNQ(I:I)
          IF(LQ .NE. ' ') THEN
            NNAME=NNAME+1
            IF(LQ .GE. 'a' .AND. LQ .LE. 'z') LQ=CHAR(ICHAR(LQ)-32)
            FNAME(NNAME:NNAME)=LQ
          ENDIF
        ENDDO
      ENDIF
      IF(JFILE .NE. 0) THEN
        OPEN(3,FILE=FNAME,STATUS='NEW',ERR=111)
        GO TO 112
111       WRITE(*,2002) ' ! '//FNAME(:NNAME)//' exists.  Overwrite?'
          READ(*,2000) LQ
          IF(LQ .NE. 'Y' .AND. LQ .NE. 'y') GO TO 11
          OPEN(3,FILE=FNAME,STATUS='UNKNOWN')
112     CONTINUE
        WRITE(3,2000) ' '
        JFILE=0
      ENDIF
      IF(QT .LE. 0.) QT=50.
      IF(QL .LE. 0.) QL=50.
      IF(QC(1) .LE. 0.) THEN
        QC(1)=50.
        QC(2)=0.
      ENDIF
      ZT=QERFI(QT/100.)
      ZL=QERFI(QL/100.)
      NC=0
      DO JC=1,4
        IF(QC(JC) .LE. 0.) GO TO 121
        NC=JC
        ZC(JC)=QERFI(QC(JC)/100.)
      ENDDO
121   CONTINUE
      IF(D0 .LE. 0.) D0=DS1
      IF(D0 .LE. 0.) D0=2.
      IF(D1 .LE. D0 .OR. DS1 .LE. 0.) THEN
        ND=1
        D1=D0
      ELSE
        DS=DS1
        ND=(D1-D0)/DS+1.75
        D1=D0+(ND-1)*DS
      ENDIF
      IF(D2 .LE. D1 .OR. DS2 .LE. 0.) THEN
        NDC=0
      ELSE
        NDC=ND
        ND=(D2-D1)/DS2+0.75
        D2=D1+ND*DS2
        ND=NDC+ND
      ENDIF
      IPOL=MAX(INDEX("HV01",CHAR(KPOL))-1,0)
      IF(IPOL .GT. 1) IPOL=IPOL-2
      JQ=MAX(INDEX("RCV012",CHAR(KST1))-1,0)
      IF(JQ .GT. 2) JQ=JQ-3
      KST(1)=JQ
      JQ=MAX(INDEX("RCV012",CHAR(KST2))-1,0)
      IF(JQ .GT. 2) JQ=JQ-3
      KST(2)=JQ
      KLIM=MAX(INDEX("CMS1234567",CHAR(KLMM))-1,0)+5
      IF(KLIM .GT. 7) KLIM=KLIM-7
      IVAR=MAX(INDEX("SAMB0123",CHAR(KVAR))-1,0)
      IF(IVAR .GT. 3) IVAR=IVAR-4
      WRITE(3,2000) TITLE,' '
      WRITE(*,2000) ' ',' ',' ',' ',' *******************************'
      WRITE(*,2000) ' ',TITLE,' '
      KWX=0
      HG(1)=SH1
      HG(2)=SH2
      DH=SDH
      CALL QLRPS(FMHZ,0.,EN0,IPOL,EPS,SGM)
      CALL QLRA(KST,KLIM,IVAR)
      WRITE(3,2015) FMHZ
      WRITE(3,2016) HG
      WRITE(3,2017) HE,KST
      WRITE(3,2018) DH
      WRITE(*,2015) FMHZ
      WRITE(*,2016) HG
      WRITE(*,2017) HE,KST
      WRITE(*,2018) DH
2015  FORMAT(12X,'Frequency',F12.0,' MHz')
2016  FORMAT(6X,'Antenna heights',2F8.1,' m')
2017  FORMAT(4X,'Effective heights',2F8.1,' m (siting=',I1,',',I1,')')
2018  FORMAT(5X,'Terrain, Delta H',F12.0,' m')
      Q=GMA/GME
      WRITE(3,2021) IPOL,EPS,SGM
      WRITE(3,2022) KLIM,ENS,Q
      WRITE(3,2000) ' '
      WRITE(*,2021) IPOL,EPS,SGM
      WRITE(*,2022) KLIM,ENS,Q
      WRITE(*,2000) ' '
2021  FORMAT(3X,'Pol=',I1,', eps=',F3.0,', sgm=',F6.3,' S/m')
2022  FORMAT(3X,'Clim=',I1,', Ns=',F5.0,'N-units, K=',F6.3)
      IF(IVAR .EQ. 0) THEN
        WRITE(3,2030)
        WRITE(*,2030)
      ELSE IF(IVAR .EQ. 1) THEN
        WRITE(3,2031) QT
        WRITE(*,2031) QT
      ELSE IF(IVAR .EQ. 2) THEN
        WRITE(3,2032) QT
        WRITE(*,2032) QT
      ELSE IF(IVAR .EQ. 3) THEN
        WRITE(3,2033) QT,QL
        WRITE(*,2033) QT,QL
      ENDIF
2030  FORMAT(3X,'Single-message service')
2031  FORMAT(3X,'Accidental service:  ',F5.1,'% time availability')
2032  FORMAT(3X,'Mobile service: required reliability-',F5.1,'%')
2033  FORMAT(3X,'Broadcast service: required reliability-',
     X            F5.1,'% time,',F5.1,'% locations')
      WRITE(3,2000) ' '
      WRITE(3,2040)
      WRITE(3,2000) ' '
      WRITE(3,2041) (QC(JC),JC=1,NC)
      WRITE(3,2000) ' '
      WRITE(*,2000) ' '
      WRITE(*,2040)
      WRITE(*,2000) ' '
      WRITE(*,2041) (QC(JC),JC=1,NC)
      WRITE(*,2000) ' '
2040  FORMAT(3X,'Estimated quantiles of basic transmission loss (dB)')
2041  FORMAT(7X,'Dist',4X,'Free',6X,'with confidence'/
     X            8X,'km',5X,'space ',4F8.1)
2045  FORMAT(2X,3F9.1,3F8.1)
      D=D0
      DT=DS
      DO JD=1,ND
        LVAR=MAX(1,LVAR)
        CALL LRPROP(D*AKM)
        FS=DB*LOG(2.*WN*DIST)
        DO JC=1,NC
          XLB(JC)=FS+AVAR(ZT,ZL,ZC(JC))
        ENDDO
        WRITE(3,2045) D,FS,(XLB(JC),JC=1,NC)
        IF(JD .EQ. 1) WRITE(*,2045) D,FS,(XLB(JC),JC=1,NC)
        IF(JD .EQ. NDC) DT=DS2
        D=D+DT
      ENDDO
      WRITE(3,2000) ' '
      WRITE(*,2000) ' '
      IF(KWX .NE. 0) THEN
        IF(KWX .EQ. 1) THEN
          WRITE(3,2081)
          WRITE(*,2081)
2081    FORMAT(
     .    3X,'!Warning! Some parameters are nearly out of range...'/
     .    3X,'   Results should be used with caution')
        ELSE IF(KWX .EQ. 2) THEN
          WRITE(3,2082)
          WRITE(*,2082)
2082    FORMAT(
     .    3X,'!Note- Default parameters have been substituted'/
     .    3X,'   for impossible ones.')
        ELSE IF(KWX .EQ. 3) THEN
          WRITE(3,2083)
          WRITE(*,2083)
2083    FORMAT(
     .    3X,'!Warning! A combination of parameters is out of range.'/
     .    3x,'   Results are probably invalid.')
        ELSE
          WRITE(3,2084)
          WRITE(*,2084)
2084    FORMAT(
     .    3X,'!Warning! Some parameters are out of range...'/
     .    3X,'   Results are probably invalid.')
        ENDIF
      ENDIF
      WRITE(3,2000) ' '
      WRITE(*,2000) ' '
      WRITE(*,2002) 'exit (N)?'
      READ(*,2000) LQ
      IF(.NOT.(LQ .EQ. 'Y' .OR. LQ .EQ. 'y')) GO TO 10
C
80    CONTINUE
        END
      SUBROUTINE VCARD(NPS,KPS,LPS,IPS,*)
        CHARACTER KPS*(*),LPS(*)*8
        INTEGER IPS(*)
C
C       A "virtual card" for input.
C       NPS, the number of parameters to be input through the
C         same number of columns.
C       KPS, a string of NPS characters indicating the corresponding
C         type of variable:  A, a single character; I, an integer;
C         digit, a real, displayed with that many decimal places.
C       LPS, an array of NPS character strings of length 8;  these
C         will be displayed as column headings, short reminders for
C         what the parameters represent.
C       IPS, an array of NPS parameter values; on entry, used to display
C         and to retain as default values; on exit, the input values.
C
C       NPS can be at most 7.  IPS is generic--integer or real--the
C         calling program must be able to reference each properly.
C         The character return is as the integer value of the ASCII
C         code of its upper case version.
C
C
      CHARACTER LLS(7)*10,L1*1,LLA*8
      CHARACTER FFMT*6
      EQUIVALENCE(IPA,APA)
C
C       The following formats are machine dependent
C       For Microsoft Fortran--
3001  FORMAT(1X,7(2X,A8))
3002  FORMAT(2H |,7A10)
3003  FORMAT(7A10)
3004  FORMAT(A,1X,$)
C
C       The following are ANSI
2000  FORMAT(A)
2001  FORMAT(I6)
2002  FORMAT(I10)
2003  FORMAT(F10.0)
C
      DATA FFMT/'(F8.X)'/
C
C
100   CONTINUE
      DO 10 I=1,NPS
      IPA=IPS(I)
      L1=KPS(I:I)
      IF(L1 .EQ. 'A') THEN
        LLA=' '
        LLA(5:5)=CHAR(IPA)
      ELSE IF(L1 .EQ. 'I') THEN
        WRITE(LLA,2001) IPA
      ELSE
        FFMT(5:5)=L1
        WRITE(LLA,FFMT) APA
      ENDIF
      LLS(I)=LLA//' |'
10    CONTINUE
      WRITE(*,3001) (LPS(I),I=1,NPS)
      WRITE(*,3002) (LLS(I),I=1,NPS)
      READ (*,3003) (LLS(I),I=1,NPS)
      DO 11 I=1,NPS
      IF(LLS(I) .NE. ' ') THEN
        L1=KPS(I:I)
        IF(L1 .EQ. 'A') THEN
          J=0
101       J=J+1
          L1=LLS(I)(J:J)
          IF(L1 .EQ. ' ') GO TO 101
          IPA=ICHAR(L1)
          IF(IPA .GE. 96) IPA=IPA-32
        ELSE IF(L1 .EQ. 'I') THEN
          READ(LLS(I),2002,ERR=13) IPA
        ELSE
          READ(LLS(I),2003,ERR=13) APA
        ENDIF
        IPS(I)=IPA
      ENDIF
11    CONTINUE
      GO TO 80
C
13    CONTINUE
      IF(I .EQ. 1) THEN
        WRITE(*,3004) '   finished?'
        READ(*,2000) L1
        IF(L1 .EQ. 'N' .OR. L1 .EQ. 'n') GO TO 100
          GO TO 81
      ENDIF
      WRITE(*,2000) '   no, no, ...'
      GO TO 100
C
81    CONTINUE
      RETURN 1
80    CONTINUE
      RETURN
        END

ext: for

max_stars_repo_name: DreamsFuture/itm

max_issues_repo_name: DreamsFuture/itm

max_forks_repo_name: DreamsFuture/itm


!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of forward_compute in forward (tangent) mode:
!   variations   of useful results: *d *e *f *concold *g *temp
!                *w *headold *x *sdata *head *dbc_data *bcperiod
!                *tempold *propunit *presold *conc *pres *a *b
!                *c doutput
!   with respect to varying inputs: *d *e *f *concold *g *temp
!                *w *headold *x *sdata *head *dbc_data *bcperiod
!                *tempold *propunit *presold *conc *pres *a *b
!                *c dinput
!   RW status of diff variables: *d:in-out *e:in-out *f:in-out
!                *concold:in-out *g:in-out *temp:in-out *w:in-out
!                *headold:in-out *x:in-out *sdata:in-out *head:in-out
!                *dbc_data:in-out *bcperiod:in-out *tempold:in-out
!                *propunit:in-out *tsal:(loc) *presold:in-out *conc:in-out
!                *pres:in-out *a:in-out *b:in-out *c:in-out doutput:out
!                dinput:in (global)thetaf:(loc) (global)thetat:(loc)
!   Plus diff mem management of: d:in e:in f:in concold:in g:in
!                temp:in w:in headold:in x:in sdata:in head:in
!                dbc_data:in bcperiod:in tempold:in propunit:in
!                tsal:in presold:in conc:in pres:in simtime:in
!                a:in b:in c:in
! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!>    @brief dependency wrapper (parameter -> simulation) useful for INVERSE and SIMUL
!>    @param[in] dinput parameter vector for initial setup (dependency vector)
!>    @param[out] doutput simulated/computed data/result values
!>    @param[in] simtime_run start time of the simulation
!>    @param[in] simtime_end finish time of the simulation
!>    @param[in] iter_out inverse iteration, SM realisation
!>    @param[in] iseed 0: FW simulation, 1 .. <mpara>: AD seeding index
!>    @param[in] ismpl local sample index
SUBROUTINE g_FORWARD_COMPUTE(dinput, g_dinput, doutput, g_doutput&
& , simtime_run, simtime_end, iter_out, iseed, ismpl)
  USE ARRAYS

  USE g_ARRAYS

  USE MOD_DATA
  USE MOD_GENRL
  IMPLICIT NONE
  INTEGER :: ismpl
  INTEGER :: i
!       input vector
  DOUBLE PRECISION :: dinput(mpara)
  DOUBLE PRECISION :: g_dinput(mpara)
!       output vector
  DOUBLE PRECISION :: doutput(ndata)
  DOUBLE PRECISION :: g_doutput(ndata)
  DOUBLE PRECISION :: trun, tend
!       parameter declaration
  INTEGER :: iter_out, iseed
  DOUBLE PRECISION :: simtime_run, simtime_end
  EXTERNAL DUMMY
!       benchmark time meassurment stuff
#ifdef BENCH
#endif
!
#ifdef  BENCH
  CALL SYS_CPUTIME(trun)
#endif
!
!       dependency setup of the parameters
  DO i=1,mpara
    CALL g_SET_OPTIP(i, dinput(i), g_dinput(i), ismpl)
  END DO
!       run forward iteration
  CALL g_FORWARD_ITER(simtime_run, simtime_end, iter_out, iseed, ismpl&
&                )
!       copy result values
  g_doutput = 0.D0
  CALL g_DCOPY(ndata, sdata(1, ismpl), g_sdata(1, ismpl), 1, doutput&
&          , g_doutput, 1)
!
!
#ifdef  BENCH
  CALL SYS_CPUTIME(tend)
  WRITE(*, '(1A,I7,1A,F9.2,1A,2F20.2)') '  [T] : component =', iseed, &
& ', time =', tend - trun, ' sec, [run] [end]:', trun, tend
#endif
!
  RETURN
END SUBROUTINE g_FORWARD_COMPUTE


text:
!
!  Written by Leandro Martínez, 2009-2011.
!  Copyright (c) 2009-2018, Leandro Martínez, Jose Mario Martinez,
!  Ernesto G. Birgin.
!
! Compute gradient relative to atom-to-atom distances
!

subroutine gparc(icart,firstjcart)

  use sizes
  use compute_data
  implicit none

  ! SCALAR ARGUMENTS
  integer :: icart,firstjcart

  ! LOCAL SCALARS
  integer :: jcart
  double precision :: datom, dtemp, xdiff, tol, &
                      short_tol, short_tol_scale

  jcart = firstjcart
  do while ( jcart .ne. 0 )
    !
    ! Cycle if this type is not to be computed
    !
    if ( .not. comptype(ibtype(jcart))) then
      jcart = latomnext(jcart)
      cycle
    end if
    !
    ! Cycle if the atoms are from the same molecule
    !
    if ( ibmol(icart) == ibmol(jcart) .and. &
         ibtype(icart) == ibtype(jcart) ) then
      jcart = latomnext(jcart)
      cycle
    end if
    !
    ! Cycle if both atoms are from fixed molecules
    !
    if ( fixedatom(icart) .and. fixedatom(jcart) ) then
      jcart = latomnext(jcart)
      cycle
    end if
    !
    ! Otherwise, compute distance and evaluate function for this pair
    !
    tol = (radius(icart)+radius(jcart))**2
    datom = (xcart(icart, 1)-xcart(jcart, 1))**2 + &
            (xcart(icart, 2)-xcart(jcart, 2))**2 + &
            (xcart(icart, 3)-xcart(jcart, 3))**2
    if( datom < tol ) then
      dtemp = fscale(icart)*fscale(jcart) * 4.d0 * (datom - tol)
      xdiff = dtemp*(xcart(icart,1) - xcart(jcart,1))
      gxcar(icart,1)= gxcar(icart,1) + xdiff
      gxcar(jcart,1)= gxcar(jcart,1) - xdiff
      xdiff = dtemp*(xcart(icart,2) - xcart(jcart,2))
      gxcar(icart,2)= gxcar(icart,2) + xdiff
      gxcar(jcart,2)= gxcar(jcart,2) - xdiff
      xdiff = dtemp*(xcart(icart,3) - xcart(jcart,3))
      gxcar(icart,3)= gxcar(icart,3) + xdiff
      gxcar(jcart,3)= gxcar(jcart,3) - xdiff
      if ( use_short_radius(icart) .or. use_short_radius(jcart) ) then
        short_tol = ( short_radius(icart) + short_radius(jcart) )**2
        if ( datom < short_tol ) then
          short_tol_scale = dsqrt(short_radius_scale(icart)*short_radius_scale(jcart))
          short_tol_scale = short_tol_scale*( tol**2 / short_tol**2 )
          dtemp = fscale(icart)*fscale(jcart) * 4.d0 * short_tol_scale*(datom - short_tol)
          xdiff = dtemp*(xcart(icart,1) - xcart(jcart,1))
          gxcar(icart,1)= gxcar(icart,1) + xdiff
          gxcar(jcart,1)= gxcar(jcart,1) - xdiff
          xdiff = dtemp*(xcart(icart,2) - xcart(jcart,2))
          gxcar(icart,2)= gxcar(icart,2) + xdiff
          gxcar(jcart,2)= gxcar(jcart,2) - xdiff
          xdiff = dtemp*(xcart(icart,3) - xcart(jcart,3))
          gxcar(icart,3)= gxcar(icart,3) + xdiff
          gxcar(jcart,3)= gxcar(jcart,3) - xdiff
        end if
      end if
    end if
    jcart = latomnext(jcart)
  end do
  return
end subroutine gparc

max_stars_repo_name: donshen/packmol

max_issues_repo_name: donshen/packmol

max_forks_repo_name: donshen/packmol