text:
!==========================================================================
elemental function gsw_ct_maxdensity (sa, p)
!==========================================================================
!
!  Calculates the Conservative Temperature of maximum density of seawater.
!  This function returns the Conservative temperature at which the density
!  of seawater is a maximum, at given Absolute Salinity, SA, and sea
!  pressure, p (in dbar).
!
!  SA =  Absolute Salinity                                         [ g/kg ]
!  p  =  sea pressure                                              [ dbar ]
!        ( i.e. absolute pressure - 10.1325 dbar )
!
!  CT_maxdensity  =  Conservative Temperature at which            [ deg C ]
!                    the density of seawater is a maximum for
!                    given Absolute Salinity and pressure.
!--------------------------------------------------------------------------

use gsw_mod_toolbox, only : gsw_alpha

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: sa, p

real (r8) :: gsw_ct_maxdensity

integer :: number_of_iterations
real (r8) :: alpha, ct, ct_mean, ct_old, dalpha_dct

real (r8), parameter :: dct = 0.001_r8

ct = 3.978_r8 - 0.22072_r8*sa         ! the initial guess of ct.

dalpha_dct = 1.1e-5_r8                ! the initial guess for dalpha_dct.

do number_of_iterations = 1, 3
    ct_old = ct
    alpha = gsw_alpha(sa,ct_old,p)
    ct = ct_old - alpha/dalpha_dct
    ct_mean = 0.5_r8*(ct + ct_old)
    dalpha_dct = (gsw_alpha(sa,ct_mean+dct,p) &
                  - gsw_alpha(sa,ct_mean-dct,p))/(dct + dct)
    ct = ct_old - alpha/dalpha_dct
end do

! After three iterations of this modified Newton-Raphson (McDougall and
! Wotherspoon, 2012) iteration, the error in CT_maxdensity is typically no
! larger than 1x10^-15 degress C.

gsw_ct_maxdensity = ct

return
end function

!--------------------------------------------------------------------------

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


subroutine restart

      use name0,   only: czero, &
                         cp9,   &
                         cp96
      use scalars, only: irst,  &
                         neqs
      use xstuff,  only: xc,    &
                         xcdot, &
                         xstore
      use time,    only: delt

      implicit none

      integer :: l

      if (irst.eq.1) then
        ! copy current state vector into backup
        call dcopy(neqs, xc, 1, xstore, 1)
      else
        ! zero velocity, restore state from backup
        do l = 1, neqs
          xcdot(l) = czero
          xc(l)    = xstore(l)
        enddo

        ! reduce time step;
        ! will arrive here with irst.eq.2 or irst.eq.3
        ! irst.eq.2 --> (irst-2) = 0, (3-irst) = 1
        !           --> delt *= 0.9
        ! irst.eq.3 --> (irst-2) = 3, (3-irst) = 0
        !           --> delt *= 0.96
        delt = delt*(cp96*(irst-2) + cp9*(3-irst))

        ! reset to normal state
        irst = 1

      endif

      return
end


text:
!##############################################################################
!# ****************************************************************************
!# <Name> transport_meshadaptation </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains all routines which are required to perform
!# mesh adaptation for scalar conservation laws in arbitrary spatial
!# dimensions.
!#
!# The following routines are available:
!#
!# 1.) transp_adaptTriangulation
!#     -> Performs h-adaptation for the given triangulation
!#
!#
!# The following auxiliary routines are available:
!#
!# 1.) transp_hadaptCallback1d
!#     -> Performs application specific tasks in the adaptation algorithm in 1D
!#
!# 2.) transp_hadaptCallback2d
!#     -> Performs application specific tasks in the adaptation algorithm in 2D
!#
!# 3.) transp_hadaptCallback3d
!#     -> Performs application specific tasks in the adaptation algorithm in 3D
!#
!# </purpose>
!##############################################################################

module transport_meshadaptation

#include "flagship.h"

!$ use omp_lib
  use basicgeometry
  use collection
  use flagship_callback
  use fsystem
  use genoutput
  use hadaptaux
  use hadaptivity
  use linearsystemscalar
  use linearsystemblock
  use triangulation

  ! Modules from transport model
  use transport_callback1d
  use transport_callback2d
  use transport_callback3d

  implicit none

  private

  public :: transp_adaptTriangulation

contains

  !*****************************************************************************

!<subroutine>

  subroutine transp_adaptTriangulation(rhadapt, rtriangulationSrc,&
      rindicator, rcollection, rtriangulationDest)

!<description>
    ! This subroutine performs h-adaptation for the given triangulation
!</description>

!<inputoutput>
    ! adaptation structure
    type(t_hadapt), intent(inout) :: rhadapt

    ! source triangulation structure
    type(t_triangulation), intent(inout), target :: rtriangulationSrc

    ! element-wise indicator
    type(t_vectorScalar), intent(inout) :: rindicator

    ! collection structure
    type(t_collection), intent(inout) :: rcollection
!</inputoutput>

!<output>
    ! OPTIONAL: destination triangulation structure
    ! If it is not given, the source triangulation is updated
    type(t_triangulation), intent(out), optional, target :: rtriangulationDest
!</output>
!</subroutine>

    ! local variables
    type(t_triangulation), pointer :: p_rtriangulation

    ! Check if adaptation structure has been prepared
    if (rhadapt%iSpec .eq. HADAPT_HAS_PARAMETERS) then

      ! Initialise adaptation structure from triangulation
      call hadapt_initFromTriangulation(rhadapt, rtriangulationSrc)

    else

      ! Refresh adaptation structure
      call hadapt_refreshAdaptation(rhadapt, rtriangulationSrc)

    end if

    ! How many spatial dimensions do we have?
    select case(rtriangulationSrc%ndim)
    case (NDIM1D)
      call transp_hadaptCallback1D(&
          HADAPT_OPR_INITCALLBACK, rcollection)
      call hadapt_performAdaptation(rhadapt, rindicator,&
          rcollection, transp_hadaptCallback1D)

    case (NDIM2D)
      call transp_hadaptCallback2D(&
          HADAPT_OPR_INITCALLBACK, rcollection)
      call hadapt_performAdaptation(rhadapt, rindicator,&
          rcollection, transp_hadaptCallback2D)

    case (NDIM3D)
      call transp_hadaptCallback3D(&
          HADAPT_OPR_INITCALLBACK, rcollection)
      call hadapt_performAdaptation(rhadapt, rindicator,&
          rcollection, transp_hadaptCallback3D)
    end select

    ! Do we have a destination triangulation or should the source
    ! triangulation structure be updated?
    if (present(rtriangulationDest)) then
      p_rtriangulation => rtriangulationDest
    else
      p_rtriangulation => rtriangulationSrc
    end if

    ! Generate raw mesh from adaptation structure
    call hadapt_generateRawMesh(rhadapt, p_rtriangulation)

  end subroutine transp_adaptTriangulation

  !*****************************************************************************

!<subroutine>

  subroutine transp_hadaptCallback1d(iOperation, rcollection)

!<description>
    ! This callback function is used to perform postprocessing tasks
    ! such as insertion/removal of elements and or vertices in the
    ! grid adaptivity procedure in 1D.
!</description>

!<input>
    ! Identifier for the grid modification operation
    integer, intent(in) :: iOperation
!</input>

!<inputoutput>
    ! A collection structure to provide additional
    ! information to the coefficient routine.
    ! This subroutine assumes the following data:
    !   rvectorQuickAccess1: solution vector
    !   IquickAccess(1):     NEQ or ivt
    !   IquickAccess(2:3):   ivt1,ivt2
    type(t_collection), intent(inout) :: rcollection
!</inputoutput>
!</subroutine>

    ! local variables
    type(t_vectorBlock), pointer, save :: rsolution
    real(DP), dimension(:), pointer, save :: p_Dsolution

    ! What operation should be performed?
    select case(iOperation)

    case(HADAPT_OPR_INITCALLBACK)
      ! Retrieve solution vector from colletion and set pointer
      rsolution => rcollection%p_rvectorQuickAccess1
      call lsysbl_getbase_double(rsolution, p_Dsolution)

      ! Call the general callback function
      call flagship_hadaptCallback1d(iOperation, rcollection)

    case(HADAPT_OPR_DONECALLBACK)
      ! Nullify solution vector
      nullify(rsolution, p_Dsolution)

      ! Call the general callback function
      call flagship_hadaptCallback1d(iOperation, rcollection)

    case(HADAPT_OPR_ADJUSTVERTEXDIM)
      ! Resize solution vector
      if (rsolution%NEQ .ne. rcollection%IquickAccess(1)) then
        call lsysbl_resizeVectorBlock(rsolution,&
            rcollection%IquickAccess(1), .false.)
        call lsysbl_getbase_double(rsolution, p_Dsolution)
      end if

    case(HADAPT_OPR_INSERTVERTEXEDGE)
      ! Insert vertex into solution vector
      if (rsolution%NEQ .lt. rcollection%IquickAccess(1)) then
        call lsysbl_resizeVectorBlock(rsolution,&
            rcollection%IquickAccess(1), .false.)
        call lsysbl_getbase_double(rsolution, p_Dsolution)
      end if
      p_Dsolution(rcollection%IquickAccess(1)) = &
          0.5_DP*(p_Dsolution(rcollection%IquickAccess(2))+&
                  p_Dsolution(rcollection%IquickAccess(3)))

      ! Call the general callback function
      call flagship_hadaptCallback1d(iOperation, rcollection)

    case(HADAPT_OPR_REMOVEVERTEX)
      ! Remove vertex from solution
      if (rcollection%IquickAccess(2) .ne. 0) then
        p_Dsolution(rcollection%IquickAccess(1)) = &
            p_Dsolution(rcollection%IquickAccess(2))
      else
        p_Dsolution(rcollection%IquickAccess(1)) = 0.0_DP
      end if

      ! Call the general callback function
      call flagship_hadaptCallback1d(iOperation, rcollection)

    case default
      ! Call the general callback function
      call flagship_hadaptCallback1d(iOperation, rcollection)

    end select

  end subroutine transp_hadaptCallback1d

  !*****************************************************************************

!<subroutine>

  subroutine transp_hadaptCallback2d(iOperation, rcollection)

!<description>
    ! This callback function is used to perform postprocessing tasks
    ! such as insertion/removal of elements and or vertices in the
    ! grid adaptivity procedure in 2D.
!</description>

!<input>
    ! Identifier for the grid modification operation
    integer, intent(in) :: iOperation
!</input>

!<inputoutput>
    ! A collection structure to provide additional
    ! information to the coefficient routine.
    ! This subroutine assumes the following data:
    !   rvectorQuickAccess1: solution vector
    !   IquickAccess(1):     NEQ or ivt
    !   IquickAccess(2:5):   ivt1,...,ivt5
    type(t_collection), intent(inout) :: rcollection
!</inputoutput>
!</subroutine>

    ! local variables
    type(t_vectorBlock), pointer, save :: rsolution
    real(DP), dimension(:), pointer, save :: p_Dsolution

    ! What operation should be performed
    select case(iOperation)

    case(HADAPT_OPR_INITCALLBACK)
      ! Retrieve solution vector from collection
      rsolution => rcollection%p_rvectorQuickAccess1
      call lsysbl_getbase_double(rsolution, p_Dsolution)

      ! Call the general callback function
      call flagship_hadaptCallback2d(iOperation, rcollection)

    case(HADAPT_OPR_DONECALLBACK)
      ! Nullify solution vector
      nullify(rsolution, p_Dsolution)

      ! Call the general callback function
      call flagship_hadaptCallback2d(iOperation, rcollection)

    case(HADAPT_OPR_ADJUSTVERTEXDIM)
      ! Resize solution vector
      if (rsolution%NEQ .ne. rcollection%IquickAccess(1)) then
        call lsysbl_resizeVectorBlock(rsolution,&
            rcollection%IquickAccess(1), .false.)
        call lsysbl_getbase_double(rsolution, p_Dsolution)
      end if

    case(HADAPT_OPR_INSERTVERTEXEDGE)
      ! Insert vertex into solution vector
      if (rsolution%NEQ .lt. rcollection%IquickAccess(1)) then
        call lsysbl_resizeVectorBlock(rsolution,&
            rcollection%IquickAccess(1), .false.)
        call lsysbl_getbase_double(rsolution, p_Dsolution)
      end if
      p_Dsolution(rcollection%IquickAccess(1)) =&
          0.5_DP*(p_Dsolution(rcollection%IquickAccess(2))+&
                  p_Dsolution(rcollection%IquickAccess(3)))

      ! Call the general callback function
      call flagship_hadaptCallback2d(iOperation, rcollection)

    case(HADAPT_OPR_INSERTVERTEXCENTR)
      ! Insert vertex into solution vector
      if (rsolution%NEQ .lt. rcollection%IquickAccess(1)) then
        call lsysbl_resizeVectorBlock(rsolution,&
            rcollection%IquickAccess(1), .false.)
        call lsysbl_getbase_double(rsolution, p_Dsolution)
      end if
      p_Dsolution(rcollection%IquickAccess(1)) =&
          0.25_DP*(p_Dsolution(rcollection%IquickAccess(2))+&
                   p_Dsolution(rcollection%IquickAccess(3))+&
                   p_Dsolution(rcollection%IquickAccess(4))+&
                   p_Dsolution(rcollection%IquickAccess(5)))

      ! Call the general callback function
      call flagship_hadaptCallback2d(iOperation, rcollection)

    case(HADAPT_OPR_REMOVEVERTEX)
      ! Remove vertex from solution
      if (rcollection%IquickAccess(2) .ne. 0) then
        p_Dsolution(rcollection%IquickAccess(1)) =&
            p_Dsolution(rcollection%IquickAccess(2))
      else
        p_Dsolution(rcollection%IquickAccess(1)) = 0.0_DP
      end if

      ! Call the general callback function
      call flagship_hadaptCallback2d(iOperation, rcollection)

    case default
      ! Call the general callback function
      call flagship_hadaptCallback2d(iOperation, rcollection)

    end select

  end subroutine transp_hadaptCallback2d

  !*****************************************************************************

!<subroutine>

  subroutine transp_hadaptCallback3d(iOperation, rcollection)

!<description>
    ! This callback function is used to perform postprocessing tasks
    ! such as insertion/removal of elements and or vertices in the
    ! grid adaptivity procedure in 3D.
!</description>

!<input>
    ! Identifier for the grid modification operation
    integer, intent(in) :: iOperation
!</input>

!<inputoutput>
    ! A collection structure to provide additional
    ! information to the coefficient routine.
    ! This subroutine assumes the following data:
    !   rvectorQuickAccess1: solution vector
    !   IquickAccess(1):     NEQ or ivt
    !   IquickAccess(2:5):   ivt1,...,ivt5
    type(t_collection), intent(inout) :: rcollection
!</inputoutput>
!</subroutine>

    ! local variables
    type(t_vectorBlock), pointer, save :: rsolution
    real(DP), dimension(:), pointer, save :: p_Dsolution

    ! What operation should be performed
    select case(iOperation)

    case(HADAPT_OPR_INITCALLBACK)
      ! Retrieve solution vector from colletion and set pointer
      rsolution => rcollection%p_rvectorQuickAccess1
      call lsysbl_getbase_double(rsolution, p_Dsolution)

      ! Call the general callback function
      call flagship_hadaptCallback3d(iOperation, rcollection)

    case(HADAPT_OPR_DONECALLBACK)
      ! Nullify solution vector
      nullify(rsolution, p_Dsolution)

      ! Call the general callback function
      call flagship_hadaptCallback3d(iOperation, rcollection)

    case(HADAPT_OPR_ADJUSTVERTEXDIM)
      ! Resize solution vector
      if (rsolution%NEQ .ne. rcollection%IquickAccess(1)) then
        call lsysbl_resizeVectorBlock(rsolution,&
            rcollection%IquickAccess(1), .false.)
        call lsysbl_getbase_double(rsolution, p_Dsolution)
      end if

    case(HADAPT_OPR_INSERTVERTEXEDGE)
      ! Insert vertex into solution vector
      if (rsolution%NEQ .lt. rcollection%IquickAccess(1)) then
        call lsysbl_resizeVectorBlock(rsolution,&
            rcollection%IquickAccess(1), .false.)
        call lsysbl_getbase_double(rsolution, p_Dsolution)
      end if
      p_Dsolution(rcollection%IquickAccess(1)) =&
          0.5_DP*(p_Dsolution(rcollection%IquickAccess(2))+&
                  p_Dsolution(rcollection%IquickAccess(3)))

      ! Call the general callback function
      call flagship_hadaptCallback3d(iOperation, rcollection)

    case(HADAPT_OPR_INSERTVERTEXCENTR)
      ! Insert vertex into solution vector
      if (rsolution%NEQ .lt. rcollection%IquickAccess(1)) then
        call lsysbl_resizeVectorBlock(rsolution,&
            rcollection%IquickAccess(1), .false.)
        call lsysbl_getbase_double(rsolution, p_Dsolution)
      end if
      p_Dsolution(rcollection%IquickAccess(1)) =&
          0.25_DP*(p_Dsolution(rcollection%IquickAccess(2))+&
                   p_Dsolution(rcollection%IquickAccess(3))+&
                   p_Dsolution(rcollection%IquickAccess(4))+&
                   p_Dsolution(rcollection%IquickAccess(5)))

      ! Call the general callback function
      call flagship_hadaptCallback3d(iOperation, rcollection)

    case(HADAPT_OPR_REMOVEVERTEX)
      ! Remove vertex from solution
      if (rcollection%IquickAccess(2) .ne. 0) then
        p_Dsolution(rcollection%IquickAccess(1)) =&
            p_Dsolution(rcollection%IquickAccess(2))
      else
        p_Dsolution(rcollection%IquickAccess(1)) = 0.0_DP
      end if

      ! Call the general callback function
      call flagship_hadaptCallback3d(iOperation, rcollection)

    case default
      ! Call the general callback function
      call flagship_hadaptCallback3d(iOperation, rcollection)
    end select

  end subroutine transp_hadaptCallback3d

end module transport_meshadaptation

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
! void c_action_intrinsic_type_spec(Token_t*, Token_t*, int, ofp_bool): Assertion `rank == 1' failed.
program p
 character(kind = 1) aK
  character(len = 2)  aL2
  character a, aL3 * 3
  character bL*2, b
  character (len = 3, kind = 1) aLK3 ! kind is lost in translation
  character(kind = 1, len = 3) aKL3 ! kind is lost in translation
  character(kind = 1) aKL4 * 4
  character(kind = 1) akL5 * (5)
  character(kind = 1, len = 7) akL6 * (6)
  character(len = 1, kind = 1) aLK
end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!  Reynolds Temp
!    NOTE:  data need to be filtered
!   output land, missing =0.

      real SST(360,180)
      real ls(360,180)
      character*80 fin,fout

!  land-sea map
      open (9,file='lstags.onedeg.dat',form='unformatted',status='old',
     *access='direct',recl=360*180*4)
      read(9,rec=1)LS
!     write(6,*)ls

!  output file
      read(5,'(a)')fout
      write(6,'(a)')fout
      open (20,file=fout,form='unformatted',status='unknown',
     *  access='direct',recl=360*180*4)
      IREC=0
      IMO=0

    1 read(5,'(a)')fin
      if(fin.eq.'END')stop
!     write(6,'(a)')fin

      open (10,file=fin,form='unformatted',status='old')

c
c  Loop for each month
  200 READ(10,END=100)
     *    IYRST,IMST,IDST,IYREND,IMEND,IDEND,NDAYS,INDEX
!     write(6,*)
!    *    IYRST,IMST,IDST,IYREND,IMEND,IDEND,NDAYS,INDEX
      READ (10) SST
      close(10)
      DO  I=1,360
        DO  J=1,180
	  if(LS(i,j).eq.0.)sst(i,j)=0.
      enddo
      enddo

      IMO =IMO + 1
      IREC=IREC+1
      write(20,rec=IREC)sst

c Print date info and SST at one location
      PRINT 7,IMO,
     * IYRST,IMST,IDST,IYREND,IMEND,IDEND,SST(70+180,80),IREC
    7  FORMAT ('Imonth =',I3,3X,'DATES =',3I4,' - ',3I4,3X,
     * 'SST (110.5W,10.5S) =',F6.2,' IREC =',i5)
       GO TO 1
c Print date info and SST at one location for last week
  100 PRINT 7,IMO,
     * IYRST,IMST,IDST,IYREND,IMEND,IDEND,SST(70+180,80),IREC

      STOP
      END

ext: f


! -----------------------------------------------------------------------------
! This file was automatically created by SARAH version SARAHVERSION
! SARAH References: arXiv:0806.0538, 0909.2863, 1002.0840, 1207.0906, 1309.7223
! (c) Florian Staub, 2013
! ------------------------------------------------------------------------------
! File created at 13:23 on 7.1.2017
! ----------------------------------------------------------------------

Module TreeLevelMasses_SM_HC

Use Control
Use Mathematics
Use MathematicsQP
Use Settings

Real(dp) :: RXiWp = 1._dp
Real(dp) :: RXiZ = 1._dp

Contains

Subroutine TreeMasses_SM(MAh,MAh2,MFd,MFd2,MFe,MFe2,MFu,MFu2,Mhh,Mhh2,MHp,               &
& MHp2,MVWp,MVWp2,MVZ,MVZ2,TW,ZDR,ZER,ZUR,ZDL,ZEL,ZUL,ZW,ZZ,v,g1,g2,g3,Lam,              &
& Yu,Yd,Ye,Mu,GenerationMixing,kont)

Implicit None

Real(dp),Intent(in) :: g1,g2,g3

Complex(dp),Intent(in) :: Lam,Yu(3,3),Yd(3,3),Ye(3,3)
Complex(dp),Intent(inout) :: Mu

Real(dp),Intent(out) :: MAh,MAh2,MFd(3),MFd2(3),MFe(3),MFe2(3),MFu(3),MFu2(3),Mhh,Mhh2,MHp,MHp2,              &
& MVWp,MVWp2,MVZ,MVZ2,TW,ZZ(2,2)

Complex(dp),Intent(out) :: ZDR(3,3),ZER(3,3),ZUR(3,3),ZDL(3,3),ZEL(3,3),ZUL(3,3),ZW(2,2)

Real(dp),Intent(in) :: v

Logical, Intent(in) :: GenerationMixing
Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4,j1,j2,j3,kontSave
Iname = Iname + 1
NameOfUnit(Iname) = 'TreeMassesSM'

kont = 0

Mu = 0.5_dp*Lam*v**2

Call CalculateVPVZ_SM(g1,g2,v,ZZ,MVZ,MVZ2,TW,kont)

Call CalculateVWp_SM(g2,v,ZW,MVWp,MVWp2,kont)

MHp2 = -1._dp*(Mu) + (v**2*Lam)/2._dp + (g2**2*v**2*RXiWp)/4._dp
  If (MHp2.ne.MHp2) Then
      Write(*,*) 'NaN appearing in MHp2'
      Call TerminateProgram
    End If
If (MHp2.Ge.0._dp) Then
 MHp = sqrt(MHp2)
Else
    MHp= 1._dp
    MHp2= 1._dp
End if

! -------------------------------
MAh2 = (-4._dp*(Mu) + 2*v**2*Lam + v**2*RXiZ*(g2*Cos(TW) + g1*Sin(TW))**2)/4._dp
  If (MAh2.ne.MAh2) Then
      Write(*,*) 'NaN appearing in MAh2'
      Call TerminateProgram
    End If
  If (MAh2.Ge.0._dp) Then
MAh = sqrt(MAh2)
  Else
        MAh= 1._dp
        MAh2= 1._dp
!    SignOfMassChanged = .True.
 End if

! -------------------------------

Mhh2 = -1._dp*(Mu) + (3*v**2*Lam)/2._dp

  If (Mhh2.ne.Mhh2) Then
      Write(*,*) 'NaN appearing in Mhh2'
      Call TerminateProgram
    End If
  If (Mhh2.Ge.0._dp) Then
Mhh = sqrt(Mhh2)
  Else
        Mhh= 1._dp
        Mhh2= 1._dp
!    SignOfMassChanged = .True.
 End if

! -------------------------------
Call CalculateMFd_SM(Yd,v,ZDL,ZDR,MFd,kont)

MFd2 = MFd**2
Call CalculateMFu_SM(Yu,v,ZUL,ZUR,MFu,kont)

MFu2 = MFu**2
Call CalculateMFe_SM(Ye,v,ZEL,ZER,MFe,kont)

MFe2 = MFe**2

 ! --------------------------------
! Setting Goldstone masses
! --------------------------------

MAh=MVZ*sqrt(RXiZ)
MAh2=MVZ2*RXiZ
MHp=MVWp*sqrt(RXiWp)
MHp2=MVWp2*RXiWp
TW = ACos(Abs(ZZ(1,1)))
Iname = Iname - 1

End Subroutine  TreeMasses_SM
!
!
! Subroutine RunningFermionMasses(MFeIN,MFe2IN,MFdIN,MFd2IN,MFuIN,MFu2IN,               &
! & v,g1,g2,g3,Lam,Yu,Yd,Ye,Mu,kont)
!
! Implicit None
!
! Integer, Intent(inout) :: kont
! Real(dp),Intent(in) :: g1,g2,g3
!
! Complex(dp),Intent(in) :: Lam,Yu(3,3),Yd(3,3),Ye(3,3),Mu
!
! Real(dp),Intent(in) :: v
!
! Real(dp),Intent(inout) :: MFeIN(3),MFe2IN(3),MFdIN(3),MFd2IN(3),MFuIN(3),MFu2IN(3)
!
! Real(dp) :: MFe(3),MFe2(3),MFd(3),MFd2(3),MFu(3),MFu2(3)
!
! Call CalculateMFe(Ye,v,ZEL,ZER,MFe,kont)
!
! MFe2 = MFe**2
! MFeIN(1:2) = MFe(1:2)
! MFe2IN(1:2) = MFe2(1:2)
! Call CalculateMFd(Yd,v,ZDL,ZDR,MFd,kont)
!
! MFd2 = MFd**2
! MFdIN(1:2) = MFd(1:2)
! MFd2IN(1:2) = MFd2(1:2)
! Call CalculateMFu(Yu,v,ZUL,ZUR,MFu,kont)
!
! MFu2 = MFu**2
! MFuIN(1:2) = MFu(1:2)
! MFu2IN(1:2) = MFu2(1:2)
! End Subroutine RunningFermionMasses
!

Subroutine TreeMassesEffPot_SM(MAh,MAh2,MFd,MFd2,MFe,MFe2,MFu,MFu2,Mhh,Mhh2,             &
& MHp,MHp2,MVWp,MVWp2,MVZ,MVZ2,TW,ZDR,ZER,ZUR,ZDL,ZEL,ZUL,ZW,ZZ,v,g1,g2,g3,              &
& Lam,Yu,Yd,Ye,Mu,GenerationMixing,kont)

Implicit None

Real(dp),Intent(in) :: g1,g2,g3

Complex(dp),Intent(in) :: Lam,Yu(3,3),Yd(3,3),Ye(3,3),Mu

Real(dp),Intent(out) :: MAh,MAh2,MFd(3),MFd2(3),MFe(3),MFe2(3),MFu(3),MFu2(3),Mhh,Mhh2,MHp,MHp2,              &
& MVWp,MVWp2,MVZ,MVZ2,TW,ZZ(2,2)

Complex(dp),Intent(out) :: ZDR(3,3),ZER(3,3),ZUR(3,3),ZDL(3,3),ZEL(3,3),ZUL(3,3),ZW(2,2)

Real(dp),Intent(in) :: v

Logical, Intent(in) :: GenerationMixing
Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4,j1,j2,j3,kontSave
Iname = Iname + 1
NameOfUnit(Iname) = 'TreeMassesSM'

kont = 0
Call CalculateVPVZEffPot_SM(g1,g2,v,ZZ,MVZ,MVZ2,TW,kont)

Call CalculateVWpEffPot_SM(g2,v,ZW,MVWp,MVWp2,kont)

MHp2 = -1._dp*(Mu) + (v**2*Lam)/2._dp + (g2**2*v**2*RXiWp)/4._dp
  If (MHp2.ne.MHp2) Then
      Write(*,*) 'NaN appearing in MHp2'
      Call TerminateProgram
    End If
  If (MHp2.Ge.0._dp) Then
MHp = sqrt(MHp2)
  Else
 End if

! -------------------------------
MAh2 = (-4._dp*(Mu) + 2*v**2*Lam + v**2*RXiZ*(g2*Cos(TW) + g1*Sin(TW))**2)/4._dp
  If (MAh2.ne.MAh2) Then
      Write(*,*) 'NaN appearing in MAh2'
      Call TerminateProgram
    End If
  If (MAh2.Ge.0._dp) Then
MAh = sqrt(MAh2)
  Else
 End if

! -------------------------------
Mhh2 = -1._dp*(Mu) + (3*v**2*Lam)/2._dp
  If (Mhh2.ne.Mhh2) Then
      Write(*,*) 'NaN appearing in Mhh2'
      Call TerminateProgram
    End If
  If (Mhh2.Ge.0._dp) Then
Mhh = sqrt(Mhh2)
  Else
 End if

!  Write(*,*) "Mhh2", Mhh2, MAh2, MHp2, MVZ2, MVWp2

! -------------------------------
Call CalculateMFdEffPot_SM(Yd,v,ZDL,ZDR,MFd,kont)

MFd2 = MFd**2
Call CalculateMFuEffPot_SM(Yu,v,ZUL,ZUR,MFu,kont)

MFu2 = MFu**2
Call CalculateMFeEffPot_SM(Ye,v,ZEL,ZER,MFe,kont)

MFe2 = MFe**2
!
!
!  If ((HighScaleModel.Eq."LOW").and.(.not.SUSYrunningFromMZ)) Then
!  If (SignOfMassChanged) Then
!  If (.Not.IgnoreNegativeMasses) Then
!   Write(*,*) " Stopping calculation because of negative mass squared."
!   Call TerminateProgram
!  Else
!   SignOfMassChanged= .False.
!   kont=0
!  End If
! End If
! If (SignOfMuChanged) Then
!  If (.Not.IgnoreMuSignFlip) Then
!   Write(*,*) " Stopping calculation because of negative mass squared in tadpoles."
!   Call TerminateProgram
!  Else
!   SignOfMuChanged= .False.
!   kont=0
!  End If
! End If
! End if
 Iname = Iname - 1

End Subroutine  TreeMassesEffPot_SM

Subroutine CalculateMFd_SM(Yd,v,ZDL,ZDR,MFd,kont)

Real(dp),Intent(in) :: v

Complex(dp),Intent(in) :: Yd(3,3)

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MFd(3)
 Complex(dp), Intent(out) :: ZDL(3,3), ZDR(3,3)

 Complex(dp) :: mat(3,3)=0._dp, mat2(3,3)=0._dp, phaseM

Complex(dp) :: ZDL2(3,3), ZDR2(3,3)

 Real(dp) :: ZDL1(3,3), ZDR1(3,3), test(2), MFd2(3)

 Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateMFd'

MFd = 0._dp
ZDL = 0._dp
ZDR = 0._dp
mat(1,1) = 0._dp
mat(1,1) = mat(1,1)+(v*Yd(1,1))/sqrt(2._dp)
mat(1,2) = 0._dp
mat(1,2) = mat(1,2)+(v*Yd(2,1))/sqrt(2._dp)
mat(1,3) = 0._dp
mat(1,3) = mat(1,3)+(v*Yd(3,1))/sqrt(2._dp)
mat(2,1) = 0._dp
mat(2,1) = mat(2,1)+(v*Yd(1,2))/sqrt(2._dp)
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)+(v*Yd(2,2))/sqrt(2._dp)
mat(2,3) = 0._dp
mat(2,3) = mat(2,3)+(v*Yd(3,2))/sqrt(2._dp)
mat(3,1) = 0._dp
mat(3,1) = mat(3,1)+(v*Yd(1,3))/sqrt(2._dp)
mat(3,2) = 0._dp
mat(3,2) = mat(3,2)+(v*Yd(2,3))/sqrt(2._dp)
mat(3,3) = 0._dp
mat(3,3) = mat(3,3)+(v*Yd(3,3))/sqrt(2._dp)

mat2 = Matmul(Transpose(Conjg(mat)),mat)
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem(Real(mat2,dp),MFd2,ZDR1,ierr,test)
ZDR2 = ZDR1
Else
Call EigenSystem(mat2,MFd2,ZDR2,ierr,test)
 End If

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(mat,Transpose(Conjg(mat)))
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem (Real(mat2,dp),MFd2,ZDL1,ierr,test)

ZDL2 = ZDL1
Else
Call EigenSystem(mat2,MFd2,ZDL2,ierr,test)

End If
ZDL2 = Conjg(ZDL2)

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(Matmul( Conjg(ZDL2),mat),Transpose( Conjg(ZDR2)))
Do i1=1,3
If (Abs(mat2(i1,i1)).gt.0._dp) Then
phaseM = mat2(i1,i1) / Abs(mat2(i1,i1))
ZDR2(i1,:) = phaseM *ZDR2(i1,:)
 End if
End Do

Do i1=1,3
If (Abs(ZDR2(i1,i1)).gt.0._dp) Then
phaseM = ZDR2(i1,i1) / Abs(ZDR2(i1,i1))
ZDR2(i1,:) = Conjg(phaseM) *ZDR2(i1,:)
 ZDL2(i1,:) = phaseM *ZDL2(i1,:)
 End if
  If (MFd2(i1).ne.MFd2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (Abs(MFd2(i1)).Le.MaxMassNumericalZero) MFd2(i1) = Abs(MFd2(i1))+1.E-10_dp
  If (MFd2(i1).Le.0._dp) Then
    If (ErrorLevel.Ge.0) Then
      Write(10,*) 'Warning from Subroutine '//NameOfUnit(Iname)
      Write(10,*) 'a mass squarred is negative: ',i1,MFd2(i1)
      Write(*,*) 'Warning from Subroutine '//NameOfUnit(Iname)
      Write(*,*) 'a mass squarred is negative: ',i1,MFd2(i1)
      Call TerminateProgram
    End If
     Write(ErrCan,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(ErrCan,*) 'in the calculation of the masses'
     Write(ErrCan,*) 'occurred a negative mass squared!'
     Write(ErrCan,*) i1,MFd2(i1)
     Write(*,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(*,*) 'in the calculation of the masses'
     Write(*,*) 'occurred a negative mass squared!'
     Write(*,*) i1,MFd2(i1)
  MFd2(i1) = 1._dp
!    SignOfMassChanged = .True.
! kont = -104
 End if
End Do

If (ierr.Ne.0.) Then
  Write(ErrCan,*) 'Warning from Subroutine CalculateMFd, ierr =',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

MFd = Sqrt( MFd2 )
ZDL = ZDL2
ZDR = ZDR2
Iname = Iname - 1

End Subroutine CalculateMFd_SM

Subroutine CalculateMFu_SM(Yu,v,ZUL,ZUR,MFu,kont)

Real(dp),Intent(in) :: v

Complex(dp),Intent(in) :: Yu(3,3)

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MFu(3)
 Complex(dp), Intent(out) :: ZUL(3,3), ZUR(3,3)

 Complex(dp) :: mat(3,3)=0._dp, mat2(3,3)=0._dp, phaseM

Complex(dp) :: ZUL2(3,3), ZUR2(3,3)

 Real(dp) :: ZUL1(3,3), ZUR1(3,3), test(2), MFu2(3)

 Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateMFu'

MFu = 0._dp
ZUL = 0._dp
ZUR = 0._dp
mat(1,1) = 0._dp
mat(1,1) = mat(1,1)-((v*Yu(1,1))/sqrt(2._dp))
mat(1,2) = 0._dp
mat(1,2) = mat(1,2)-((v*Yu(2,1))/sqrt(2._dp))
mat(1,3) = 0._dp
mat(1,3) = mat(1,3)-((v*Yu(3,1))/sqrt(2._dp))
mat(2,1) = 0._dp
mat(2,1) = mat(2,1)-((v*Yu(1,2))/sqrt(2._dp))
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)-((v*Yu(2,2))/sqrt(2._dp))
mat(2,3) = 0._dp
mat(2,3) = mat(2,3)-((v*Yu(3,2))/sqrt(2._dp))
mat(3,1) = 0._dp
mat(3,1) = mat(3,1)-((v*Yu(1,3))/sqrt(2._dp))
mat(3,2) = 0._dp
mat(3,2) = mat(3,2)-((v*Yu(2,3))/sqrt(2._dp))
mat(3,3) = 0._dp
mat(3,3) = mat(3,3)-((v*Yu(3,3))/sqrt(2._dp))

mat2 = Matmul(Transpose(Conjg(mat)),mat)
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem(Real(mat2,dp),MFu2,ZUR1,ierr,test)
ZUR2 = ZUR1
Else
Call EigenSystem(mat2,MFu2,ZUR2,ierr,test)
 End If

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(mat,Transpose(Conjg(mat)))
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem (Real(mat2,dp),MFu2,ZUL1,ierr,test)

ZUL2 = ZUL1
Else
Call EigenSystem(mat2,MFu2,ZUL2,ierr,test)

End If
ZUL2 = Conjg(ZUL2)

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(Matmul( Conjg(ZUL2),mat),Transpose( Conjg(ZUR2)))
Do i1=1,3
If (Abs(mat2(i1,i1)).gt.0._dp) Then
phaseM = mat2(i1,i1) / Abs(mat2(i1,i1))
ZUR2(i1,:) = phaseM *ZUR2(i1,:)
 End if
End Do

Do i1=1,3
If (Abs(ZUR2(i1,i1)).gt.0._dp) Then
phaseM = ZUR2(i1,i1) / Abs(ZUR2(i1,i1))
ZUR2(i1,:) = Conjg(phaseM) *ZUR2(i1,:)
 ZUL2(i1,:) = phaseM *ZUL2(i1,:)
 End if
  If (MFu2(i1).ne.MFu2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (Abs(MFu2(i1)).Le.MaxMassNumericalZero) MFu2(i1) = Abs(MFu2(i1))+1.E-10_dp
  If (MFu2(i1).Le.0._dp) Then
    If (ErrorLevel.Ge.0) Then
      Write(10,*) 'Warning from Subroutine '//NameOfUnit(Iname)
      Write(10,*) 'a mass squarred is negative: ',i1,MFu2(i1)
      Write(*,*) 'Warning from Subroutine '//NameOfUnit(Iname)
      Write(*,*) 'a mass squarred is negative: ',i1,MFu2(i1)
      Call TerminateProgram
    End If
     Write(ErrCan,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(ErrCan,*) 'in the calculation of the masses'
     Write(ErrCan,*) 'occurred a negative mass squared!'
     Write(ErrCan,*) i1,MFu2(i1)
     Write(*,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(*,*) 'in the calculation of the masses'
     Write(*,*) 'occurred a negative mass squared!'
     Write(*,*) i1,MFu2(i1)
  MFu2(i1) = 1._dp
!    SignOfMassChanged = .True.
! kont = -104
 End if
End Do

If (ierr.Ne.0.) Then
  Write(ErrCan,*) 'Warning from Subroutine CalculateMFu, ierr =',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

MFu = Sqrt( MFu2 )
ZUL = ZUL2
ZUR = ZUR2
Iname = Iname - 1

End Subroutine CalculateMFu_SM

Subroutine CalculateMFe_SM(Ye,v,ZEL,ZER,MFe,kont)

Real(dp),Intent(in) :: v

Complex(dp),Intent(in) :: Ye(3,3)

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MFe(3)
 Complex(dp), Intent(out) :: ZEL(3,3), ZER(3,3)

 Complex(dp) :: mat(3,3)=0._dp, mat2(3,3)=0._dp, phaseM

Complex(dp) :: ZEL2(3,3), ZER2(3,3)

 Real(dp) :: ZEL1(3,3), ZER1(3,3), test(2), MFe2(3)

 Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateMFe'

MFe = 0._dp
ZEL = 0._dp
ZER = 0._dp
mat(1,1) = 0._dp
mat(1,1) = mat(1,1)+(v*Ye(1,1))/sqrt(2._dp)
mat(1,2) = 0._dp
mat(1,2) = mat(1,2)+(v*Ye(2,1))/sqrt(2._dp)
mat(1,3) = 0._dp
mat(1,3) = mat(1,3)+(v*Ye(3,1))/sqrt(2._dp)
mat(2,1) = 0._dp
mat(2,1) = mat(2,1)+(v*Ye(1,2))/sqrt(2._dp)
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)+(v*Ye(2,2))/sqrt(2._dp)
mat(2,3) = 0._dp
mat(2,3) = mat(2,3)+(v*Ye(3,2))/sqrt(2._dp)
mat(3,1) = 0._dp
mat(3,1) = mat(3,1)+(v*Ye(1,3))/sqrt(2._dp)
mat(3,2) = 0._dp
mat(3,2) = mat(3,2)+(v*Ye(2,3))/sqrt(2._dp)
mat(3,3) = 0._dp
mat(3,3) = mat(3,3)+(v*Ye(3,3))/sqrt(2._dp)

mat2 = Matmul(Transpose(Conjg(mat)),mat)
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem(Real(mat2,dp),MFe2,ZER1,ierr,test)
ZER2 = ZER1
Else
Call EigenSystem(mat2,MFe2,ZER2,ierr,test)
 End If

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(mat,Transpose(Conjg(mat)))
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem (Real(mat2,dp),MFe2,ZEL1,ierr,test)

ZEL2 = ZEL1
Else
Call EigenSystem(mat2,MFe2,ZEL2,ierr,test)

End If
ZEL2 = Conjg(ZEL2)

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(Matmul( Conjg(ZEL2),mat),Transpose( Conjg(ZER2)))
Do i1=1,3
If (Abs(mat2(i1,i1)).gt.0._dp) Then
phaseM = mat2(i1,i1) / Abs(mat2(i1,i1))
ZER2(i1,:) = phaseM *ZER2(i1,:)
 End if
End Do

Do i1=1,3
If (Abs(ZER2(i1,i1)).gt.0._dp) Then
phaseM = ZER2(i1,i1) / Abs(ZER2(i1,i1))
ZER2(i1,:) = Conjg(phaseM) *ZER2(i1,:)
 ZEL2(i1,:) = phaseM *ZEL2(i1,:)
 End if
  If (MFe2(i1).ne.MFe2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (Abs(MFe2(i1)).Le.MaxMassNumericalZero) MFe2(i1) = Abs(MFe2(i1))+1.E-10_dp
  If (MFe2(i1).Le.0._dp) Then
    If (ErrorLevel.Ge.0) Then
      Write(10,*) 'Warning from Subroutine '//NameOfUnit(Iname)
      Write(10,*) 'a mass squarred is negative: ',i1,MFe2(i1)
      Write(*,*) 'Warning from Subroutine '//NameOfUnit(Iname)
      Write(*,*) 'a mass squarred is negative: ',i1,MFe2(i1)
      Call TerminateProgram
    End If
     Write(ErrCan,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(ErrCan,*) 'in the calculation of the masses'
     Write(ErrCan,*) 'occurred a negative mass squared!'
     Write(ErrCan,*) i1,MFe2(i1)
     Write(*,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(*,*) 'in the calculation of the masses'
     Write(*,*) 'occurred a negative mass squared!'
     Write(*,*) i1,MFe2(i1)
  MFe2(i1) = 1._dp
!    SignOfMassChanged = .True.
! kont = -104
 End if
End Do

If (ierr.Ne.0.) Then
  Write(ErrCan,*) 'Warning from Subroutine CalculateMFe, ierr =',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

MFe = Sqrt( MFe2 )
ZEL = ZEL2
ZER = ZER2
Iname = Iname - 1

End Subroutine CalculateMFe_SM

Subroutine CalculateVPVZ_SM(g1,g2,v,ZZ,MVZ,MVZ2,TW,kont)

Real(dp), Intent(in) :: g1,g2,v

Real(dp), Intent(out) :: TW

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MVZ, MVZ2
Real(dp) :: VPVZ2(2),VPVZ(2)

Real(dp), Intent(out) :: ZZ(2,2)

Real(dp) :: mat(2,2)

Real(dp) ::  test(2)

Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateVPVZ'

mat(1,1) = 0._dp
mat(1,1) = mat(1,1)+(g1**2*v**2)/4._dp
mat(1,2) = 0._dp
mat(1,2) = mat(1,2)-(g1*g2*v**2)/4._dp
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)+(g2**2*v**2)/4._dp

 Do i1=2,2
  Do i2 = 1, i1-1
  mat(i1,i2) = mat(i2,i1)
  End do
End do

Call EigenSystem(mat,VPVZ2,ZZ,ierr,test)

ZZ = Transpose(ZZ)
If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

If ((ierr.Ne.0.).And.(ErrorLevel.Ge.-1)) Then
  Write(10,*) 'Warning from Subroutine '//NameOfUnit(Iname)
  Write(10,*) 'Diagonalization failed, ierr : ',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

Do i1=1,2
  If (Abs(VPVZ2(i1)).Le.1.E-10_dp*(Maxval(VPVZ2))) VPVZ2(i1) = 1.E-10_dp
  If (VPVZ2(i1).ne.VPVZ2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (VPVZ2(i1).Ge.0._dp) Then
  VPVZ(i1) =Sqrt(VPVZ2(i1) )
  Else
    If (ErrorLevel.Ge.0) Then
      Write(10,*) 'Warning from Subroutine '//NameOfUnit(Iname)
      Write(10,*) 'a mass squarred is negative: ',i1,VPVZ2(i1)
    End If
  VPVZ(i1)= 1._dp
  VPVZ2(i1)= 1._dp
     Write(ErrCan,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(ErrCan,*) 'in the calculation of the masses'
     Write(ErrCan,*) 'occurred a negative mass squared!'
     Write(ErrCan,*) i1,VPVZ2(i1)
     Write(*,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(*,*) 'in the calculation of the masses'
     Write(*,*) 'occurred a negative mass squared!'
     Write(*,*) i1,VPVZ2(i1)
  VPVZ(i1)= 1._dp
  VPVZ2(i1) = 1._dp
!    SignOfMassChanged = .True.
! kont = -104
 End if
End Do

MVZ = VPVZ(2)
MVZ2 = VPVZ2(2)
TW = ACos(Abs(ZZ(1,1)))

 Iname = Iname - 1

End Subroutine CalculateVPVZ_SM

Subroutine CalculateVWp_SM(g2,v,ZW,MVWp,MVWp2,kont)

Real(dp), Intent(in) :: g2,v

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MVWp, MVWp2
Real(dp) :: VWp2(2),VWp(2)

Complex(dp), Intent(out) :: ZW(2,2)

Complex(dp) :: mat(2,2)

Real(dp) ::  test(2)

Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateVWp'

mat(1,1) = 0._dp
mat(1,1) = mat(1,1)+(g2**2*v**2)/4._dp
mat(1,2) = 0._dp
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)+(g2**2*v**2)/4._dp

 Do i1=2,2
  Do i2 = 1, i1-1
  mat(i1,i2) = Conjg(mat(i2,i1))
  End do
End do

Call EigenSystem(mat,VWp2,ZW,ierr,test)

ZW = Transpose(ZW)
If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

If ((ierr.Ne.0.).And.(ErrorLevel.Ge.-1)) Then
  Write(10,*) 'Warning from Subroutine '//NameOfUnit(Iname)
  Write(10,*) 'Diagonalization failed, ierr : ',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

Do i1=1,2
  If (Abs(VWp2(i1)).Le.1.E-10_dp*(Maxval(VWp2))) VWp2(i1) = 1.E-10_dp
  If (VWp2(i1).ne.VWp2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (VWp2(i1).Ge.0._dp) Then
  VWp(i1) =Sqrt(VWp2(i1) )
  Else
    If (ErrorLevel.Ge.0) Then
      Write(10,*) 'Warning from Subroutine '//NameOfUnit(Iname)
      Write(10,*) 'a mass squarred is negative: ',i1,VWp2(i1)
    End If
  VWp(i1)= 1._dp
  VWp2(i1)= 1._dp
     Write(ErrCan,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(ErrCan,*) 'in the calculation of the masses'
     Write(ErrCan,*) 'occurred a negative mass squared!'
     Write(ErrCan,*) i1,VWp2(i1)
     Write(*,*) 'Warning from routine '//NameOfUnit(Iname)
     Write(*,*) 'in the calculation of the masses'
     Write(*,*) 'occurred a negative mass squared!'
     Write(*,*) i1,VWp2(i1)
  VWp(i1)= 1._dp
  VWp2(i1) = 1._dp
!    SignOfMassChanged = .True.
! kont = -104
 End if
End Do

MVWp = VWp(1)
MVWp2 = VWp2(1)

 Iname = Iname - 1

End Subroutine CalculateVWp_SM
!
Subroutine CalculateMFdEffPot_SM(Yd,v,ZDL,ZDR,MFd,kont)

Real(dp),Intent(in) :: v

Complex(dp),Intent(in) :: Yd(3,3)

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MFd(3)
 Complex(dp), Intent(out) :: ZDL(3,3), ZDR(3,3)

 Complex(dp) :: mat(3,3)=0._dp, mat2(3,3)=0._dp, phaseM

Complex(dp) :: ZDL2(3,3), ZDR2(3,3)

 Real(dp) :: ZDL1(3,3), ZDR1(3,3), test(2), MFd2(3)

 Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateMFd'

MFd = 0._dp
ZDL = 0._dp
ZDR = 0._dp
mat(1,1) = 0._dp
mat(1,1) = mat(1,1)+(v*Yd(1,1))/sqrt(2._dp)
mat(1,2) = 0._dp
mat(1,2) = mat(1,2)+(v*Yd(2,1))/sqrt(2._dp)
mat(1,3) = 0._dp
mat(1,3) = mat(1,3)+(v*Yd(3,1))/sqrt(2._dp)
mat(2,1) = 0._dp
mat(2,1) = mat(2,1)+(v*Yd(1,2))/sqrt(2._dp)
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)+(v*Yd(2,2))/sqrt(2._dp)
mat(2,3) = 0._dp
mat(2,3) = mat(2,3)+(v*Yd(3,2))/sqrt(2._dp)
mat(3,1) = 0._dp
mat(3,1) = mat(3,1)+(v*Yd(1,3))/sqrt(2._dp)
mat(3,2) = 0._dp
mat(3,2) = mat(3,2)+(v*Yd(2,3))/sqrt(2._dp)
mat(3,3) = 0._dp
mat(3,3) = mat(3,3)+(v*Yd(3,3))/sqrt(2._dp)

mat2 = Matmul(Transpose(Conjg(mat)),mat)
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem(Real(mat2,dp),MFd2,ZDR1,ierr,test)
ZDR2 = ZDR1
Else
Call EigenSystem(mat2,MFd2,ZDR2,ierr,test)
 End If

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(mat,Transpose(Conjg(mat)))
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem (Real(mat2,dp),MFd2,ZDL1,ierr,test)

ZDL2 = ZDL1
Else
Call EigenSystem(mat2,MFd2,ZDL2,ierr,test)

End If
ZDL2 = Conjg(ZDL2)

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(Matmul( Conjg(ZDL2),mat),Transpose( Conjg(ZDR2)))
Do i1=1,3
If (Abs(mat2(i1,i1)).gt.0._dp) Then
phaseM = mat2(i1,i1) / Abs(mat2(i1,i1))
ZDR2(i1,:) = phaseM *ZDR2(i1,:)
 End if
End Do

Do i1=1,3
If (Abs(ZDR2(i1,i1)).gt.0._dp) Then
phaseM = ZDR2(i1,i1) / Abs(ZDR2(i1,i1))
ZDR2(i1,:) = Conjg(phaseM) *ZDR2(i1,:)
 ZDL2(i1,:) = phaseM *ZDL2(i1,:)
 End if
  If (MFd2(i1).ne.MFd2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (Abs(MFd2(i1)).Le.MaxMassNumericalZero) MFd2(i1) = Abs(MFd2(i1))+1.E-10_dp
  If (MFd2(i1).Le.0._dp) Then
! kont = -104
 End if
End Do

If (ierr.Ne.0.) Then
  Write(ErrCan,*) 'Warning from Subroutine CalculateMFd, ierr =',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

MFd = Sqrt( MFd2 )
ZDL = ZDL2
ZDR = ZDR2
Iname = Iname - 1

End Subroutine CalculateMFdEffPot_SM

Subroutine CalculateMFuEffPot_SM(Yu,v,ZUL,ZUR,MFu,kont)

Real(dp),Intent(in) :: v

Complex(dp),Intent(in) :: Yu(3,3)

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MFu(3)
 Complex(dp), Intent(out) :: ZUL(3,3), ZUR(3,3)

 Complex(dp) :: mat(3,3)=0._dp, mat2(3,3)=0._dp, phaseM

Complex(dp) :: ZUL2(3,3), ZUR2(3,3)

 Real(dp) :: ZUL1(3,3), ZUR1(3,3), test(2), MFu2(3)

 Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateMFu'

MFu = 0._dp
ZUL = 0._dp
ZUR = 0._dp
mat(1,1) = 0._dp
mat(1,1) = mat(1,1)-((v*Yu(1,1))/sqrt(2._dp))
mat(1,2) = 0._dp
mat(1,2) = mat(1,2)-((v*Yu(2,1))/sqrt(2._dp))
mat(1,3) = 0._dp
mat(1,3) = mat(1,3)-((v*Yu(3,1))/sqrt(2._dp))
mat(2,1) = 0._dp
mat(2,1) = mat(2,1)-((v*Yu(1,2))/sqrt(2._dp))
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)-((v*Yu(2,2))/sqrt(2._dp))
mat(2,3) = 0._dp
mat(2,3) = mat(2,3)-((v*Yu(3,2))/sqrt(2._dp))
mat(3,1) = 0._dp
mat(3,1) = mat(3,1)-((v*Yu(1,3))/sqrt(2._dp))
mat(3,2) = 0._dp
mat(3,2) = mat(3,2)-((v*Yu(2,3))/sqrt(2._dp))
mat(3,3) = 0._dp
mat(3,3) = mat(3,3)-((v*Yu(3,3))/sqrt(2._dp))

mat2 = Matmul(Transpose(Conjg(mat)),mat)
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem(Real(mat2,dp),MFu2,ZUR1,ierr,test)
ZUR2 = ZUR1
Else
Call EigenSystem(mat2,MFu2,ZUR2,ierr,test)
 End If

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(mat,Transpose(Conjg(mat)))
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem (Real(mat2,dp),MFu2,ZUL1,ierr,test)

ZUL2 = ZUL1
Else
Call EigenSystem(mat2,MFu2,ZUL2,ierr,test)

End If
ZUL2 = Conjg(ZUL2)

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(Matmul( Conjg(ZUL2),mat),Transpose( Conjg(ZUR2)))
Do i1=1,3
If (Abs(mat2(i1,i1)).gt.0._dp) Then
phaseM = mat2(i1,i1) / Abs(mat2(i1,i1))
ZUR2(i1,:) = phaseM *ZUR2(i1,:)
 End if
End Do

Do i1=1,3
If (Abs(ZUR2(i1,i1)).gt.0._dp) Then
phaseM = ZUR2(i1,i1) / Abs(ZUR2(i1,i1))
ZUR2(i1,:) = Conjg(phaseM) *ZUR2(i1,:)
 ZUL2(i1,:) = phaseM *ZUL2(i1,:)
 End if
  If (MFu2(i1).ne.MFu2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (Abs(MFu2(i1)).Le.MaxMassNumericalZero) MFu2(i1) = Abs(MFu2(i1))+1.E-10_dp
  If (MFu2(i1).Le.0._dp) Then
! kont = -104
 End if
End Do

If (ierr.Ne.0.) Then
  Write(ErrCan,*) 'Warning from Subroutine CalculateMFu, ierr =',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

MFu = Sqrt( MFu2 )
ZUL = ZUL2
ZUR = ZUR2
Iname = Iname - 1

End Subroutine CalculateMFuEffPot_SM

Subroutine CalculateMFeEffPot_SM(Ye,v,ZEL,ZER,MFe,kont)

Real(dp),Intent(in) :: v

Complex(dp),Intent(in) :: Ye(3,3)

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MFe(3)
 Complex(dp), Intent(out) :: ZEL(3,3), ZER(3,3)

 Complex(dp) :: mat(3,3)=0._dp, mat2(3,3)=0._dp, phaseM

Complex(dp) :: ZEL2(3,3), ZER2(3,3)

 Real(dp) :: ZEL1(3,3), ZER1(3,3), test(2), MFe2(3)

 Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateMFe'

MFe = 0._dp
ZEL = 0._dp
ZER = 0._dp
mat(1,1) = 0._dp
mat(1,1) = mat(1,1)+(v*Ye(1,1))/sqrt(2._dp)
mat(1,2) = 0._dp
mat(1,2) = mat(1,2)+(v*Ye(2,1))/sqrt(2._dp)
mat(1,3) = 0._dp
mat(1,3) = mat(1,3)+(v*Ye(3,1))/sqrt(2._dp)
mat(2,1) = 0._dp
mat(2,1) = mat(2,1)+(v*Ye(1,2))/sqrt(2._dp)
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)+(v*Ye(2,2))/sqrt(2._dp)
mat(2,3) = 0._dp
mat(2,3) = mat(2,3)+(v*Ye(3,2))/sqrt(2._dp)
mat(3,1) = 0._dp
mat(3,1) = mat(3,1)+(v*Ye(1,3))/sqrt(2._dp)
mat(3,2) = 0._dp
mat(3,2) = mat(3,2)+(v*Ye(2,3))/sqrt(2._dp)
mat(3,3) = 0._dp
mat(3,3) = mat(3,3)+(v*Ye(3,3))/sqrt(2._dp)

mat2 = Matmul(Transpose(Conjg(mat)),mat)
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem(Real(mat2,dp),MFe2,ZER1,ierr,test)
ZER2 = ZER1
Else
Call EigenSystem(mat2,MFe2,ZER2,ierr,test)
 End If

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(mat,Transpose(Conjg(mat)))
If (Maxval(Abs(Aimag(mat2))).Eq.0._dp) Then
Call EigenSystem (Real(mat2,dp),MFe2,ZEL1,ierr,test)

ZEL2 = ZEL1
Else
Call EigenSystem(mat2,MFe2,ZEL2,ierr,test)

End If
ZEL2 = Conjg(ZEL2)

If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

mat2 = Matmul(Matmul( Conjg(ZEL2),mat),Transpose( Conjg(ZER2)))
Do i1=1,3
If (Abs(mat2(i1,i1)).gt.0._dp) Then
phaseM = mat2(i1,i1) / Abs(mat2(i1,i1))
ZER2(i1,:) = phaseM *ZER2(i1,:)
 End if
End Do

Do i1=1,3
If (Abs(ZER2(i1,i1)).gt.0._dp) Then
phaseM = ZER2(i1,i1) / Abs(ZER2(i1,i1))
ZER2(i1,:) = Conjg(phaseM) *ZER2(i1,:)
 ZEL2(i1,:) = phaseM *ZEL2(i1,:)
 End if
  If (MFe2(i1).ne.MFe2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (Abs(MFe2(i1)).Le.MaxMassNumericalZero) MFe2(i1) = Abs(MFe2(i1))+1.E-10_dp
  If (MFe2(i1).Le.0._dp) Then
! kont = -104
 End if
End Do

If (ierr.Ne.0.) Then
  Write(ErrCan,*) 'Warning from Subroutine CalculateMFe, ierr =',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

MFe = Sqrt( MFe2 )
ZEL = ZEL2
ZER = ZER2
Iname = Iname - 1

End Subroutine CalculateMFeEffPot_SM

Subroutine CalculateVPVZEffPot_SM(g1,g2,v,ZZ,MVZ,MVZ2,TW,kont)

Real(dp), Intent(in) :: g1,g2,v

Real(dp), Intent(out) :: TW

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MVZ, MVZ2
Real(dp) :: VPVZ2(2),VPVZ(2)

Real(dp), Intent(out) :: ZZ(2,2)

Real(dp) :: mat(2,2)

Real(dp) ::  test(2)

Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateVPVZ'

mat(1,1) = 0._dp
mat(1,1) = mat(1,1)+(g1**2*vFix**2)/4._dp
mat(1,2) = 0._dp
mat(1,2) = mat(1,2)-(g1*g2*vFix**2)/4._dp
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)+(g2**2*vFix**2)/4._dp

 Do i1=2,2
  Do i2 = 1, i1-1
  mat(i1,i2) = mat(i2,i1)
  End do
End do

Call EigenSystem(mat,VPVZ2,ZZ,ierr,test)

ZZ = Transpose(ZZ)
If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

If ((ierr.Ne.0.).And.(ErrorLevel.Ge.-1)) Then
  Write(10,*) 'Warning from Subroutine '//NameOfUnit(Iname)
  Write(10,*) 'Diagonalization failed, ierr : ',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

Do i1=1,2
  If (Abs(VPVZ2(i1)).Le.1.E-10_dp*(Maxval(VPVZ2))) VPVZ2(i1) = 1.E-10_dp
  If (VPVZ2(i1).ne.VPVZ2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (VPVZ2(i1).Ge.0._dp) Then
  VPVZ(i1) =Sqrt(VPVZ2(i1) )
  Else
  VPVZ(i1)= 1._dp
  VPVZ(i1)= 1._dp
! kont = -104
 End if
End Do

MVZ = VPVZ(2)
MVZ2 = VPVZ2(2)

 Iname = Iname - 1

End Subroutine CalculateVPVZEffPot_SM

Subroutine CalculateVWpEffPot_SM(g2,v,ZW,MVWp,MVWp2,kont)

Real(dp), Intent(in) :: g2,v

Integer, Intent(inout) :: kont
Integer :: i1,i2,i3,i4, ierr
Integer :: j1,j2,j3,j4
Real(dp), Intent(out) :: MVWp, MVWp2
Real(dp) :: VWp2(2),VWp(2)

Complex(dp), Intent(out) :: ZW(2,2)

Complex(dp) :: mat(2,2)

Real(dp) ::  test(2)

Iname = Iname + 1
NameOfUnit(Iname) = 'CalculateVWp'

mat(1,1) = 0._dp
mat(1,1) = mat(1,1)+(g2**2*vFix**2)/4._dp
mat(1,2) = 0._dp
mat(2,2) = 0._dp
mat(2,2) = mat(2,2)+(g2**2*vFix**2)/4._dp

 Do i1=2,2
  Do i2 = 1, i1-1
  mat(i1,i2) = Conjg(mat(i2,i1))
  End do
End do

Call EigenSystem(mat,VWp2,ZW,ierr,test)

ZW = Transpose(ZW)
If ((ierr.Eq.-8).Or.(ierr.Eq.-9)) Then
  Write(ErrCan,*) "Possible numerical problem in "//NameOfUnit(Iname)
  If (ErrorLevel.Eq.2) Then
  Write(*,*) "Possible numerical problem in "//NameOfUnit(Iname)
    Call TerminateProgram
  End If
  ierr = 0
End If

If ((ierr.Ne.0.).And.(ErrorLevel.Ge.-1)) Then
  Write(10,*) 'Warning from Subroutine '//NameOfUnit(Iname)
  Write(10,*) 'Diagonalization failed, ierr : ',ierr
  kont = ierr
  Iname = Iname - 1
  Return
End If

Do i1=1,2
  If (Abs(VWp2(i1)).Le.1.E-10_dp*(Maxval(VWp2))) VWp2(i1) = 1.E-10_dp
  If (VWp2(i1).ne.VWp2(i1)) Then
      Write(*,*) 'NaN appearing in '//NameOfUnit(Iname)
      Call TerminateProgram
    End If
  If (VWp2(i1).Ge.0._dp) Then
  VWp(i1) =Sqrt(VWp2(i1) )
  Else
  VWp(i1)= 1._dp
  VWp(i1)= 1._dp
! kont = -104
 End if
End Do

MVWp = VWp(1)
MVWp2 = VWp2(1)

 Iname = Iname - 1

End Subroutine CalculateVWpEffPot_SM
!
! Subroutine TreeMassesSM(MFd,MFd2,MFe,MFe2,MFu,MFu2,MVWp,MVWp2,MVZ,MVZ2,               &
! & TW,ZDR,ZER,ZUR,ZDL,ZEL,ZUL,ZW,ZZ,v,g1,g2,g3,Lam,Yu,Yd,Ye,Mu,GenerationMixing,kont)
!
! Implicit None
!
! Real(dp),Intent(in) :: g1,g2,g3
!
! Complex(dp),Intent(in) :: Lam,Yu(3,3),Yd(3,3),Ye(3,3),Mu
!
! Real(dp),Intent(out) :: MFd(3),MFd2(3),MFe(3),MFe2(3),MFu(3),MFu2(3),MVWp,MVWp2,MVZ,MVZ2,TW,ZZ(2,2)
!
! Complex(dp),Intent(out) :: ZDR(3,3),ZER(3,3),ZUR(3,3),ZDL(3,3),ZEL(3,3),ZUL(3,3),ZW(2,2)
!
! Real(dp),Intent(in) :: v
!
! Logical, Intent(in) :: GenerationMixing
! Integer, Intent(inout) :: kont
! Integer :: i1,i2,i3,i4,j1,j2,j3,kontSave
! Iname = Iname + 1
! NameOfUnit(Iname) = 'TreeMassesSM'
!
! kont = 0
! Call CalculateMFd(Yd,v,ZDL,ZDR,MFd,kont)
!
! MFd2 = MFd**2
! Call CalculateMFu(Yu,v,ZUL,ZUR,MFu,kont)
!
! MFu2 = MFu**2
! Call CalculateMFe(Ye,v,ZEL,ZER,MFe,kont)
!
! MFe2 = MFe**2
!
!
!  Call CalculateVPVZ(g1,g2,v,ZZ,MVZ,MVZ2,TW,kont)
!
! Call CalculateVWp(g2,v,ZW,MVWp,MVWp2,kont)
!
! Iname = Iname - 1
!
! End Subroutine  TreeMassesSM
!
!
! Subroutine SortGoldstones(MAh,MAh2,MFd,MFd2,MFe,MFe2,MFu,MFu2,Mhh,Mhh2,               &
! & MHp,MHp2,MVWp,MVWp2,MVZ,MVZ2,TW,ZDR,ZER,ZUR,ZDL,ZEL,ZUL,ZW,ZZ,kont)
!
! Real(dp),Intent(inout) :: MAh,MAh2,MFd(3),MFd2(3),MFe(3),MFe2(3),MFu(3),MFu2(3),Mhh,Mhh2,MHp,MHp2,              &
! & MVWp,MVWp2,MVZ,MVZ2,TW,ZZ(2,2)
!
! Complex(dp),Intent(inout) :: ZDR(3,3),ZER(3,3),ZUR(3,3),ZDL(3,3),ZEL(3,3),ZUL(3,3),ZW(2,2)
!
! Integer, Intent(inout) :: kont
! Integer :: i1, i2, pos
!
!
! End subroutine SortGoldstones

End Module TreeLevelMasses_SM_HC


text:
The Lucille Youmans and Dixie Benches Bench Memorials honors the memory of Lucille Youmans and her dog, Dixie.  Those who sit here can enjoy a shady, pleasant view of the Aspen Greenbelt Dogs dog Dog Parks run.

Cool

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


    Program zgebal_example

!     ZGEBAL Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use blas_interfaces, Only: dznrm2
      Use lapack_example_aux, Only: nagf_sort_realvec_rank, &
        nagf_file_print_matrix_complex_gen_comp, &
        nagf_sort_cmplxvec_rank_rearrange
      Use lapack_interfaces, Only: zgebak, zgebal, zgehrd, zhseqr, ztrevc, &
        zunghr
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Complex (Kind=dp) :: scal
      Integer :: i, ifail, ihi, ilo, info, k, lda, ldh, ldvl, ldvr, lwork, m, &
        n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), h(:, :), tau(:), vl(:, :), &
        vr(:, :), vr_row(:), w(:), work(:)
      Real (Kind=dp), Allocatable :: rwork(:), scale(:)
      Integer, Allocatable :: irank(:)
      Logical :: select(1)
      Character (1) :: clabs(1), rlabs(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs, aimag, conjg, maxloc, real
!     .. Executable Statements ..
      Write (nout, *) 'ZGEBAL Example Program Results'
      Flush (nout)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n
      ldvl = 1
      lda = n
      ldh = n
      ldvr = n
      lwork = 64*n
      Allocate (a(lda,n), h(ldh,n), tau(n), vl(ldvl,1), vr(ldvr,n), w(n), &
        work(lwork), rwork(n), scale(n), irank(n))

!     Read A from data file
      Read (nin, *)(a(i,1:n), i=1, n)

!     Balance A
      Call zgebal('Both', n, a, lda, ilo, ihi, scale, info)

!     Reduce A to upper Hessenberg form H = (Q**H)*A*Q
      Call zgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)

!     Copy A to H and VR
      h(1:n, 1:n) = a(1:n, 1:n)
      vr(1:n, 1:n) = a(1:n, 1:n)

!     Form Q explicitly, storing the result in VR
      Call zunghr(n, 1, n, vr, ldvr, tau, work, lwork, info)

!     Calculate the eigenvalues and Schur factorization of A
      Call zhseqr('Schur form', 'Vectors', n, ilo, ihi, h, ldh, w, vr, ldvr, &
        work, lwork, info)

      Write (nout, *)
      If (info>0) Then
        Write (nout, *) 'Failure to converge.'
      Else
!       Sort eigenvalues into descending absolute value
        rwork(1:n) = abs(w(1:n))
!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call nagf_sort_realvec_rank(rwork, 1, n, 'Descending', irank, ifail)
!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call nagf_sort_cmplxvec_rank_rearrange(w, 1, n, irank, ifail)

        Write (nout, *) 'Eigenvalues'
        Write (nout, 100)(' (', real(w(i)), ',', aimag(w(i)), ')', i=1, n)
        Flush (nout)

!       Calculate the eigenvectors of A, storing the result in VR
        Call ztrevc('Right', 'Backtransform', select, n, h, ldh, vl, ldvl, vr, &
          ldvr, n, m, work, rwork, info)

        Call zgebak('Both', 'Right', n, ilo, ihi, scale, m, vr, ldvr, info)

!       Reorder eigenvectors using irank
        Allocate (vr_row(n))
        Do i = 1, n
          vr_row(1:n) = vr(i, 1:n)
!         ifail: behaviour on error exit
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
          ifail = 0
          Call nagf_sort_cmplxvec_rank_rearrange(vr_row, 1, n, irank, ifail)
          vr(i, 1:n) = vr_row(1:n)
        End Do
        Deallocate (vr_row)

!       Print eigenvectors

        Write (nout, *)
        Flush (nout)

!       Normalize the eigenvectors, largest element real
        Do i = 1, m
          rwork(1:n) = abs(vr(1:n,i))
          k = maxloc(rwork(1:n), 1)
          scal = conjg(vr(k,i))/abs(vr(k,i))/dznrm2(n, vr(1,i), 1)
          vr(1:n, i) = vr(1:n, i)*scal
        End Do

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call nagf_file_print_matrix_complex_gen_comp('General', ' ', n, m, vr, &
          ldvr, 'Bracketed', 'F7.4', 'Contents of array VR', 'Integer', rlabs, &
          'Integer', clabs, 80, 0, ifail)

      End If

100   Format ((3X,4(A,F7.4,A,F7.4,A,:)))
    End Program


text:
************************************************************************
* Read parameters for time discretization
*
* This routine reads the parameters for the time discretization
* and stores them in the corresponding COMMON blockdiscretization.
*
* In:
*   MDATA  - Unit number to use for reading process
*   MSHOW  - Level of output in the initialization phase
*   CFNAME - Name of the file
************************************************************************

      SUBROUTINE RDTIDC (MDATA,MSHOW,CFNAME)

      IMPLICIT NONE

      INCLUDE 'cout.inc'
      INCLUDE 'cfileout.inc'
      INCLUDE 'cerr.inc'

      INCLUDE 'ctimediscr.inc'

C     parameters

      CHARACTER CFNAME*(*)
      INTEGER MSHOW,MDATA

C     local variables

      INTEGER I,IFMTS
      CHARACTER CSTR*(255)

C     Open DAT file

      IFMTS = 1
      CALL OF0 (MDATA,CFNAME,IFMTS)

C     Ignore the first 6 lines

      DO I=1,6
        READ (MDATA,*)
      END DO

C     Check version number

      READ (MDATA,*) I
      IF (I.NE.100) THEN
        WRITE (*,'(A)')
     *        'Error in reading time-discretization parameters'
        WRITE (*,'(A)') 'Version number incorrect'
        STOP
      END IF

      WRITE (CSTR,'(A)') ' Time discretization parameters'
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)
      WRITE (CSTR,'(A)') '--------------------------------'
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     Ignore separator

      READ (MDATA,*)

C     Stationary or nonstationary prolem

      READ (MDATA,*) ISTAT
      IF ((ISTAT.LT.0).OR.(ISTAT.GT.1)) ISTAT=0
      WRITE (CSTR,1000)
     *      'Time dependency               : ISTAT  = ', ISTAT
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     Separator

      READ (MDATA,*)
      WRITE (MTERM,9001)

C     Mesh reconstruction/adaption in time

      READ (MDATA,*) IDMESH
      WRITE (CSTR,'(A,I10)')
     *      'Mesh reconstruction/adaption:   IDMESH = ', IDMESH
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) IMSCRT
      WRITE (CSTR,'(A,I10)')
     *      'Mesh adaption criterion:        IMSCRT = ', IMSCRT
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) IMSREP
      WRITE (CSTR,'(A,I10)')
     *      'Adaption repetitions:           IMSREP = ', IMSREP
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) IMSSTP
      WRITE (CSTR,'(A,I10)')
     *      'Timesteps between mesh adapt.:  IMSSTP = ', IMSSTP
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) DMSMTM
      WRITE (CSTR,1001)
     *      'Max. Time for grid adaption   : DMSMTM = ', DMSMTM
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     Separator

      READ (MDATA,*)
      WRITE (MTERM,9001)

C     Stopping criteria on time-dependent solver

      READ (MDATA,*) NITNS
      WRITE (CSTR,'(A,I10)')
     *      'Max. Number of time steps     : NITNS  = ', NITNS
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) TIMEMX
      WRITE (CSTR,1001)
     *      'Max. Time                     : TIMEMX = ', TIMEMX
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) TIMENS
      WRITE (CSTR,1001)
     *      'Absolute start time           : TIMENS = ', TIMENS
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) TSTEP
      WRITE (CSTR,1001)
     *      '(Initial) time step           : TSTEP  = ', TSTEP
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) IFRSTP
      WRITE (CSTR,1000)
     *      'Time stepping scheme          : IFRSTP = ', IFRSTP
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) THETA
      WRITE (CSTR,1001)
     *      'THETA for one-step scheme     : THETA  = ', THETA
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) IEPSAD
      IF ((IEPSAD.LT.0).OR.(IEPSAD.GT.8)) IEPSAD=1
      WRITE (CSTR,1000)
     *      'Acceptance criterion          : IEPSAD = ', IEPSAD
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) EPSNS
      WRITE (CSTR,1001)
     *      'Limit for time derivative     : EPSNS  = ', EPSNS
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     Separator

      READ (MDATA,*)
      WRITE (MTERM,9001)

C     Parameters for the adaptive time stepping

      READ (MDATA,*) IADTIM
      IF (ABS(IADTIM).GT.3) IADTIM=0
      WRITE (CSTR,1000)
     *      'Type of adaptivity            : IADTIM = ', IADTIM
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) IREPIT
      IREPIT = ABS(IREPIT)
      WRITE (CSTR,1000)
     *      'Max.numbers of repetitions    : IREPIT = ', IREPIT
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) DTMIN
      WRITE (CSTR,1001)
     *      'Min. Timestep                 : DTMIN  = ', DTMIN
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) DTMAX
      WRITE (CSTR,1001)
     *      'Max. Timestep                 : DTMAX  = ', DTMAX
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) DTFACT
      DTFACT=ABS(DTFACT)
      IF (DTFACT.LT.1D0) DTFACT=1D0
      WRITE (CSTR,1001)
     *      'Max. Timestep change          : DTFACT = ', DTFACT
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) EPSADU
      WRITE (CSTR,1001)
     *      'EPS for not acceptance        : EPSADU = ', EPSADU
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     Separator

      READ (MDATA,*)
      WRITE (MTERM,9001)

C     Start procedure

      READ (MDATA,*) IADIN
      IF ((IADIN.LT.0).OR.(IADIN.GT.2)) IADIN=0
      WRITE (CSTR,1000)
     *      'Type of start procedure       : IADIN  = ', IADIN
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) TIMEIN
      WRITE (CSTR,1001)
     *      'Time length of start procedure: TIMEIN = ', TIMEIN
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) EPSADI
      WRITE (CSTR,1001)
     *      'EPS for accep. in start proc. : EPSADI = ', EPSADI
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) EPSADL
      WRITE (CSTR,1001)
     *      'EPS for accept. after start pr: EPSADL = ', EPSADL
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     That's it.

      WRITE (MTERM,*)

      CLOSE (MDATA)

1000  FORMAT (A,I4)
1001  FORMAT (A,E16.8)
1002  FORMAT (A,2I4)
1003  FORMAT (A,A)

9000  FORMAT(79('-'))
9001  FORMAT(60('-'))

      END


text:
! mpiexec -np 4 bin/Debug/WGFM_test2.exe
program mainWGF_solve
    use meshread_mod
    use linalg_mod
    use WGFM_matrices
    use spherePEC_mod
    use dipole_functions
    use planewave_sources
    use dipole_mod
    use currents_mod
    use mpi
    implicit none

    real(8), parameter :: pi = 4d0 * atan(1d0)
    complex(8), parameter :: IU = (0d0, 1d0)

    ! Variables
    type(Mesh) :: msh

    complex(8), allocatable :: ZMat(:, :)  ! WGFM matrix

    complex(8), allocatable :: Esrc(:, :), Hsrc(:, :)
    complex(8), allocatable :: MJsrc(:)
    complex(8), allocatable :: Uvec(:), Vvec(:), UVvec(:)

    ! MPI variables
    integer :: ierr, N_procs, id

    ! Parameters
    real(8), parameter :: ep0 = 1d0 / (35950207149.4727056d0 * PI)  ! vacuum permittivity [F/m]
    real(8), parameter :: mu0 = PI * 4d-7                           ! vacuum permeability [H/m]
    real(8), parameter :: k0 = 0.5d0          ! free-space wavenumber [1/m]

    real(8), parameter :: ep1_r = 1d0         ! upper half-space relative permittivity
    real(8), parameter :: ep2_r = 7d0         ! lower half-space relative permittivity
    real(8), parameter :: ep1 = ep0*ep1_r     ! upper half-space permittivity [F/m]
    real(8), parameter :: mu1 = mu0           ! upper half-space permeability [H/m]
    real(8), parameter :: ep2 = ep0*ep2_r     ! lower half-space permittivity [F/m]
    real(8), parameter :: mu2 = mu0           ! lower half-space permeability [H/m]

    real(8), parameter :: w = k0 / sqrt(ep0 * mu0)   ! angular frequency [rad/s]
    real(8), parameter :: k1 = w * sqrt(ep1 * mu1)  ! upper half-space wavenumber [1/m]
    real(8), parameter :: k2 = w * sqrt(ep2 * mu2)  ! lower half-space wavenumber [1/m]

    character(len=100), parameter :: file_msh = 'meshes/plane_mesh/plane_mesh_a40_h14.msh'
    real(8), parameter :: win_radius = 20d0   ! Window radius

    ! Dipole/planewave parameters
    real(8), parameter :: pol(3) = [0d0, 1d0, 0d0]    ! dipole polarization vector
    real(8), parameter :: src(3) = [0d0, 0d0, 5d0]    ! location of dipole source

    !************************************************************************************
    ! MPI init
    call mpi_init(ierr)
    call mpi_comm_size(MPI_COMM_WORLD, N_procs, ierr)
    call mpi_comm_rank(MPI_COMM_WORLD, id, ierr)
    print *, "Total processes:", N_procs, "Process id:", id

    ! Set window type and size
    call set_circ_window(win_radius, c_in=0.7d0)

    ! Load mesh
    call load_gmsh(file_msh, msh)

    ! Compute WGFM matrix
    if (id == 0) print *, "WGFM"
    call genWGFMMat(msh, w, k1, k2, ep1, ep2, mu1, mu2, ZMat)
    if (id == 0) print *, "done"

    ! Compute source fields
    if (id == 0) then
        ! Source fields
        call ElectricDipoleOnGmsh(pol, src, k1, w, mu1, msh, Esrc, Hsrc)

        ! WGFM RHS
        call genWGFMRHS(msh, Esrc, Hsrc, MJsrc)

        ! Solve
        print *, "solving"
        allocate(Uvec(msh%nbEdg))
        allocate(Vvec(msh%nbEdg))
        allocate(UVvec(2*msh%nbEdg))
        call linsolve(UVvec, ZMat, MJsrc, 2*msh%nbEdg)
        Uvec = UVvec(1 : msh%nbEdg)
        Vvec = UVvec(msh%nbEdg+1 : 2*msh%nbEdg)
        deallocate(MJsrc, UVvec, ZMat)
        print *, "done"

        ! Save currents
        call save_msh_currents(file_msh, msh, k1, k2, w,  &
                               ep1_r, ep2_r, win_radius, Uvec, Vvec)
    end if

    call MPI_Finalize(ierr)
end program

max_stars_repo_name: Riarrieta/WindowedGreenFunctionMethod

max_issues_repo_name: Riarrieta/WindowedGreenFunctionMethod

max_forks_repo_name: Riarrieta/WindowedGreenFunctionMethod


text:
The Fitnessgarage Sports Academy offers massage massage therapy and performance training focusing on core strength for youth and adult sports athletes. Their core training program emphasizes healthful lifestyle, strength training, sports conditioning, flexibility and agility. Core strength focuses on strengthening the abdominals, back, hamstrings and hip flexors to provide better flexibility, strength and long, lean musculature. Unlike most gyms, there is no initiation fee, nor will you find state of the art specialized exercise equipment  the focus of this facility is on fitness training rather than body building. Private and group lessons are available for youth and adults.

The Fitnessgarage hosts a variety of classes and services to tailor to your needs (see their http://fitnessgarage.us/schedule.html classes schedule).  Their program incorporates what we call CORETraining, a process that includes Performance Training, and Tutoring.

Most of the classes are $10 per 1hour session, and you can get a punch card for 5 classes for $45.  Classes tend to be fairly small with quite a bit of personal attention.  Hideshi is easily able to work around different needsfor example, Hideshi has combined shoulder with lower back and core exercises during the Friday morning lower back class, and has indicated that if someone who needs to work on the lower back drops by one of the shoulder classes, he can easily incorporate it into that hour.

The owner, Hideshi Okamoto, and many of his staff have a track & field background. Okamoto is a coach for the Japanese Olympic track and field team, a Japanese track sports commentator, a board member of the National Association of Speed and Explosion and a certified massage therapist.

 PHOTO ALBUMS

http://s854.photobucket.com/albums/ab105/fitnessgarage/Lacrosse%20Nov%203%202010/ LACROSSE TRAINING 2010 NOV. 3
https://picasaweb.google.com/103506024750631640639/SpeedAndAgilityCamp2011Jun?authuser0&featdirectlink Speed and Agility Camp 2011 Jun.
https://picasaweb.google.com/103506024750631640639/SpeedAgilityCamp2011Aug1st5th?authuser0&featdirectlink Speed and Agility Camp 2011 Aug.1st5th
https://picasaweb.google.com/103506024750631640639/SpeedAndAgilityCamp2011August812?authuser0&featdirectlink Speed and Agility Camp 2011 Aug.8th12th

Therapeutic massage
Open Hour
Mon. to Fri.
10:00am4:00pm
Need appointment Call or Email us

Massage Fee
Half hour session $30
1 hour session $60

 Classes
Hyumas Kick Boxercise
Each class is $10/session
A CARDIO based balance challenge for adults!
Date Every Wed.
Time 8:159:45
Only accepting 8 people.

Shoulder Injury Prevention Class
Each class is $10/session
Wed. 7:008:00am

Lower Back Injury Prevention Class
Each class is $10/session
Fri. 8:309:30am

Low Impact Class
Each class is $10/session
Time: 11:0012:00
Date: Mon. / Wed. / Fri.

Core and Balance
These classes are designed to focus on improving ones balance.By focusing on balance, we are able to strengthen our smaller stabilizer muscle, which limit your joint movement while helping the bigger muscles (ones that do the lifting) to complete the bodys movement.Balance exercises improve ones range of motion (joint mobilization and muscle stiffness) which decreases and as we age.

Each class is $10/session
Mon. 8:309:30 am
Tue.  6:007:00 am / 8:159:45 am
Wed. 8:009:00 am
Thu.  6:007:00 am / 8:159:45 am

Saturday Morning HardCore
One of the most fascinating muscle groups in your body, YOUR CORE. Also referred to as your power house. We are firm believers that one is only as strong as his or her core. The core include your front abs (rectus abdominis), side abs (obliques), lower back, and inner thigh muscles. While all our classes do include at least ten minutes devoted to this area, it is the fastest recovering area of your body hence the hour and a half session. Please do not feel as though you need to stay the entire 90 minutes! This is more of a rotating class where you can come in anytime, perferrably before 9:30 and stay for an hour or as long as you like.

Each class is $10/session
Sat.  9:0010:30 am

Youth Class
Our youth classes focus on developing the core through a series of exercises based on balance, speed, agility,coordination and Strength. Open to all ages(717), all abilities, and all sports.
Call about a group pricing plan. Office (530) 2971736, Cell (530) 2193807

Each class is $20/session
or $60/month (once a week)
or $100/month(twice a week)

Youth Class Punch Card is available now!
5times for $90

Mon. 4:005:00 pm / 5:006:00 pm
Tue. 4:005:00 pm / 5:006:00 pm / 6:007:00 pm
Wed.4:005:00 pm
Thu. 4:005:00 pm/5:006:00 pm / 6:007:00 pm
Fri.  4:005:00 pm / 5:006:00 pm

We have more classes, (see their http://fitnessgarage.us/schedule.html classes schedule)

 Special Class

Lacrosse PreSeason Speed & Agility Training

Session16weeks($120)
9/12 thru 10/19 Mon & Wed
4pm5pm or 5pm6pm

Session24weeks($80)
10/24 thru 11/16 Mon & Wed
4pm5pm or 5pm6pm

http://fitnessgarage.us/web/lacrosse_2011.pdf More information
http://fitnessgarage.us/web/lacrosse_fall%20_2011.pdf Application Form

DHS Baseball Team Training

10/1112/1 (8weeks $160)
Tue @ 3PM & Thr @ 2pm

http://fitnessgarage.us/web/baseball2011.pdf Application Form

Adult Running Class by Hideshi Okamoto
Start @ Monday, September 12 · 6:00pm
Each class is $10

Place:Rainbow Park(behind Arts Center on F st.& 14th)
Date: Mon. 6pm/Wed. 9am

Sprint Clinic by former Olympic Coach (88, 92, 96, 00) Hideshi Okamoto
Start @ 4:30pm
for Junior and High school Kids
(@4:00pm for Elementally kids)
Call us! Limited space Max 10Kids!

  See our Gyms & Fitness page to find out about other gyms in the Davis area.

20110118 14:04:14 nbsp Hideshi has been my neighbor and friend for several years.  Several times he has fixed bad shoulder pain with his massage and manipulations.  He does excellent work.  Users/JimStewart

20110527 09:59:52 nbsp I just got back from the 1hour lower back injury prevention class, having experienced intermittent lower back trouble for years.  Hideshi was fantastic, and its a really nice setup there.  The place is back in the corner of the parking lot, not really visible from the street, but its not that hard to find.  Theres no big gym equipment there; the training is dynamic and constantly active.  If youre looking to bulk up, this isnt the place to do it.  There are no memberships.  Each class you take is $10, and if they recognize you, youre welcome to drop in at other times of the day and use the equipment.

Im somewhat out of shape because of a recent injury, but not too bad.  The exercises were fairly low intensity, but constant.  I left with my quads and hamstrings quite tired and slightly shaky, my hip flexors, lower back, and abs all feeling well worked, and my upper back, shoulders, tris, and bis all feeling pleasantly tired.  I dont know whether Ill be able to work it into my work schedule, but Id love to go back as frequently as possible.  Im also looking forward to trying out the Tuesday/Thursday core and balance class.

I should also note that Hideshi himself is awesome.  Extremely friendly and engaging and a great instructor. Users/TomGarberson

The fitness garage is nestled away in East Davis, but Hideshi is worth every penny and then some Daubert


text:
      subroutine cmo_create(cmo_name,ierror_return)
C
C
C#######################################################################
C
C      PURPOSE -
C
C         This Routine Creates a Mesh Object.
C
C      INPUT ARGUMENTS -
C
C         cmo_name - (character) Mesh Object Name.
C
C      OUTPUT ARGUMENTS -
C
C         ierror_return - Error Return Code (==0 ==> OK, <>0 ==> Error)
C
C      CHANGE HISTORY -
C
C         $Log: cmo_create.f,v $
C         Revision 2.00  2007/11/05 19:45:48  spchu
C         Import to CVS
C
CPVCS
CPVCS       Rev 1.1   Thu Jan 20 16:08:08 2000   dcg
CPVCS
CPVCS       Rev 1.24   Fri Nov 05 13:26:56 1999   dcg
CPVCS    remove dictionary dependencies
CPVCS
CPVCS       Rev 1.23   Mon Apr 14 16:40:08 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.22   Mon Nov 18 13:57:52 1996   dcg
CPVCS    define sbname='sbcmoprm' correctly
CPVCS
CPVCS       Rev 1.21   Mon Nov 18 10:38:20 1996   dcg
CPVCS    remove character literals from argument lists
CPVCS
CPVCS       Rev 1.20   09/11/95 14:43:24   het
CPVCS    Change to the storage block based CMO stuff.
CPVCS
CPVCS       Rev 1.19   08/30/95 21:07:30   het
CPVCS    Put cmo table data into the cmoatt storage block
CPVCS
CPVCS       Rev 1.18   08/29/95 12:49:10   dcg
CPVCS    set length for names to 40 characters
CPVCS
CPVCS       Rev 1.17   08/29/95 12:16:04   het
CPVCS    Add the cmoatt storage block for each CMO
CPVCS
CPVCS       Rev 1.16   08/22/95 06:51:52   het
CPVCS    Split the storage block for CMO variables.
CPVCS
CPVCS       Rev 1.15   05/23/95 06:50:26   het
CPVCS    Change dictionary so that they are CMO specific
CPVCS
CPVCS       Rev 1.14   03/15/95 15:22:34   ejl
CPVCS    Finished installing the defaults.
CPVCS
CPVCS       Rev 1.13   02/16/95 09:55:52   ejl
CPVCS    Fixed bugs, fixed hole in the Create command.
CPVCS    Added commands MODATT, LENGTH, MEMORY, & COMPRESS.
CPVCS
CPVCS       Rev 1.12   02/10/95 14:06:44   ejl
CPVCS    Fix bugs left from last update.
CPVCS
CPVCS       Rev 1.10   01/30/95 18:17:44   het
CPVCS    Add the cmo_select call to refresh cmo.h pointers
CPVCS
CPVCS       Rev 1.9   01/30/95 06:22:08   het
CPVCS    Fix several cmo errors
CPVCS
CPVCS       Rev 1.8   01/24/95 08:52:26   het
CPVCS    Add error checking to the cmo routines.
CPVCS
CPVCS
CPVCS       Rev 1.7   01/23/95 17:02:32   het
CPVCS    Correct some character problems with the table driven
CPVCS       cmo attributes.
CPVCS
CPVCS
CPVCS       Rev 1.6   01/23/95 12:38:06   het
CPVCS    Use the table driven attribute management changes.
CPVCS
CPVCS
CPVCS       Rev 1.5   01/04/95 22:01:32   llt
CPVCS    unicos changes (made by het)
CPVCS
CPVCS       Rev 1.4   01/02/95 12:36:38   het
CPVCS
CPVCS
CPVCS       Rev 1.3   01/02/95 12:26:36   het
CPVCS
CPVCS
CPVCS       Rev 1.2   12/21/94 09:50:02   het
CPVCS    Add a call to 'cmo_select()' to select the just created
CPVCS       cmo as the current_cmo.
CPVCS
CPVCS       Rev 1.1   12/11/94 17:51:48   het
CPVCS    Fixed error related to cmo_create.
CPVCS
CPVCS
CPVCS       Rev 1.0   12/09/94 22:49:32   het
CPVCS    Original version.
CPVCS
C
C#######################################################################
C
      implicit none
C
C#######################################################################
C
      include 'cmo_lg.h'
C
C#######################################################################
C
      character*(*) cmo_name
C
      integer ierror_return
C
C#######################################################################
C
C     LOCAL VARIABLE DEFINITION
C
      integer len
C
      integer ierr, icmonum
C
      character*32 isubname
C
      character*132 logmess
C
      integer icscode, iexist
C
C#######################################################################
C
      integer icharlnf
C
C#######################################################################
C
C
C.... Define the routine name.
C
      isubname='cmo_create'
C
C.... Check name against reserved names.
C
      len=icharlnf(cmo_name)
C
      if((cmo_name(1:len).eq.'-cmo-') .or.
     *   (cmo_name(1:len).eq.'-def-') .or.
     *   (cmo_name(1:len).eq.'-all-') .or.
     *   (cmo_name(1:len).eq.'-default-') .or.
     *   (cmo_name(1:len).eq.'-notset-')) then
C
C....    Name is reserved.
C
         ierror_return=-1
         write(logmess,'(a,a)')
     *      '    ERROR: Mesh Object name is reserved: ',cmo_name
         call writloga('default',0,logmess,0,ierr)
C
      endif
C
C
C.... Search table for Mesh Object.
C
      call cmo_exist(cmo_name,iexist)
C
      if(iexist.eq.0) then
C
C....    Mesh Object already exists.
C
         ierror_return=-1
         write(logmess,'(a,a)')
     *      'Mesh Object already exists: ',cmo_name
         call writloga('default',0,logmess,0,ierr)
C
      else
C
C....    Build a directory entry for the new Mesh Object.
C
         call cmo_create_dir(cmo_name,ierror_return)
C
C
C....    Set the default attributes.
C
         call cmo_set_default(cmo_name,icmonum,ierror_return)
C
C
C....    Set up the new Mesh Object.
C
         call cmo_allocate(cmo_name,ierror_return)
c
C....    Make this the Current Mesh Object.
C
         call cmo_select(cmo_name,ierror_return)
C
      endif
C
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


text:
!
! This file is released under terms of BSD license
! See LICENSE file for more information
!
! Simple program to test the loop-fusion directive
!

PROGRAM LOOP_FUSION
  CALL clawloop
END

SUBROUTINE clawloop
  INTEGER :: i, k

  !$acc parallel
  !$acc loop seq
  DO k=1,10
    !$acc loop
    !$claw loop-fusion
    DO i=1,10
      PRINT *, 'First loop body:',i
    END DO

    !$acc loop
    !$claw loop-fusion
    DO i=1,10
      PRINT *, 'Second loop body:',i
    END DO
  END DO
  !$acc end parallel

END

max_stars_repo_name: peclatj/claw-compiler

max_issues_repo_name: peclatj/claw-compiler

max_forks_repo_name: peclatj/claw-compiler