! initilize movement types
SUBROUTINE movetype_init()
  use movetype
  IMPLICIT NONE
  integer :: ierr
  write(*,*) "movetype_init:"
  ALLOCATE(movetype_name(nmovetypes),stat=ierr)
  ALLOCATE(movetype_prob(nmovetypes),stat=ierr)
  ALLOCATE(movetype_ntry(nmovetypes),stat=ierr)
  ALLOCATE(movetype_nsuccess(nmovetypes),stat=ierr)
  ALLOCATE(movetype_i_success(nmovetypes),stat=ierr)
  ALLOCATE(movetype_i_try(nmovetypes),stat=ierr)
  movetype_ntry = 0
  movetype_nsuccess = 0
  movetype_i_success = 0
  movetype_i_try = 0
END SUBROUTINE

! initilize movement types
SUBROUTINE movetype_set_ensemble(name)
  use movetype
  use coupling_pres
  use coupling_exch
  IMPLICIT NONE
  CHARACTER(LEN=5) :: name
  integer :: ind_str

  write(*,*) "movetype_set_ensemble:"
  ensemble_temp = .False.
  ensemble_pres = .False.
  ensemble_exch = .False.
  name = adjustl(name)

  nmovetypes = 0
  IF (scan(name,"T") /= 0) THEN ! if ensemble_name contain "T"
    ensemble_temp = .True.
    write(*,*) " Turn on temperature coupling."
    nmovetypes = nmovetypes + 1
  ENDIF
  If (scan(name,"P") /= 0) THEN
    ensemble_pres = .True.
    nmovetypes = nmovetypes + 1
    if (scan(name,"Z") /= 0) THEN
      semiiso = 3
      write(*,*) " Turn on semiisotropic pressure coupling on z-direction."
      write(*,*) " But, not supported yet"
      stop
    else
      semiiso = 0
      write(*,*) " Turn on isotropic pressure coupling."
    endif
  ENDIF
  IF (scan(name,"X") /= 0) THEN
    ensemble_exch = .True.
    nmovetypes = nmovetypes + 1
    !ind_str = scan(name,"X")+1
    !exch_comp = name(ind_str:ind_str)
    !write(*,*) " Turn on exchange identity coupling based on the component, ", exch_comp
    write(*,*) " Turn on exchange identity coupling"
  ENDIF
  if (nmovetypes == 0) THEN
    WRITE(*,*) " no coupling? because of your ", name, " ensemble."
    STOP
  ENDIF
  write(*,*) " total #couplings = ", nmovetypes
  return
END SUBROUTINE movetype_set_ensemble

SUBROUTINE movetype_read(fileunit)
  use movetype
  use trans
  use coupling_pres
  use coupling_exch

  IMPLICIT NONE
  CHARACTER(len=5) :: coupling, exch_opt_run
  integer :: fileunit, i, j, ierr
  double precision :: temp, tprob

  write(*,*) "movetype_read:"
  do i=1, nmovetypes
    read(fileunit,*) coupling
    write(*,*) coupling
    ! translation
    if (trim(coupling) == 'trans') THEN
      movetype_name(i) = 'trans'
      read(fileunit,*) movetype_prob(i), dlr_a, dlr_b ! probability, translation increment for components
      write(*,'(1x,A,1X,3(F10.5),1x,A)') " trans = ",movetype_prob(i), dlr_a, dlr_b, "(prob, dlr_a, dlr_b)"
      cycle
    ! pressure coupling
    else if (trim(coupling) == 'press') THEN
      movetype_name(i) = 'press'
      read(fileunit,*) movetype_prob(i), press_val, dvx, temp ! probability, pressure, volume change, temperature
      write(*,'(1x,A,1X,4(F10.5),1x,A)') " press = ", movetype_prob(i), press_val, dvx, temp, "(prob, press., volume change, temp.)"
      tinv = 1.0/temp
      cycle
    ! particle exchange coupling
    else if (trim(coupling) == 'exch') THEN
      movetype_name(i) = 'exch'
      read(fileunit,*) movetype_prob(i), xi_val, exch_ncomp, exch_opt_run ! probability, xi value, number of component
      if(exch_opt_run == 'OPT') then
        xi_val_opt = .true.
        write(*,*) " Turn on optimizing xi_val."
      else
        xi_val_opt = .false.
      endif
      write(*,'(1x,A,1X,2(F10.5),I5,1x,A)') " exch = ", movetype_prob(i), xi_val, exch_ncomp, "(prob, xi_val, exchange component)"
      if(exch_ncomp <= 1) then
        write(*,*) " exch_ncomp should be greater than one component"
        stop
      else if (xi_val >= 1.0) then
        write(*,*) " xi_val, fugacity ratio, should be less than one"
        stop
      end if
      log_xi = dlog(xi_val/(1.0D0-xi_val))
      ALLOCATE(exch_tcomp(exch_ncomp),stat=ierr)
      !ALLOCATE(exch_typs(nptot),stat=ierr)
      exch_typs = 0 ! initialize
      read(fileunit,*) (exch_tcomp(j), j=1, exch_ncomp)
      write(*,*) " component => ", (exch_tcomp(j), j=1, exch_ncomp)
      write(*,'(A,A,A,A,A,A)') " From your setting, component change xi_(", exch_tcomp(1), &
                "), attempts from species ", exch_tcomp(1), " to species ", exch_tcomp(2)
      write(*,*) ", which gives m = -1. Otherwise, m = 1"
      write(*,'(A,A,A)') " (in other words, xi_val is xi_", exch_tcomp(1),")"
      write(*,*) "you set log(xi/(1-xi)) = ", log_xi
      cycle
    else
      write(*,*) " wrong some lines"
      stop
    endif
  enddo
  tprob = 0.0D0
  do i=1, nmovetypes
    tprob = tprob + movetype_prob(i)
  enddo
  if (tprob < 1.0D0 .or. tprob > 1.0D0) then
    write(*,*) " total probability is not 1, ", tprob, ". check your input file."
    stop
  endif
  return
END SUBROUTINE movetype_read

! exit movetype
SUBROUTINE movetype_exit
  use movetype
  use coupling_exch
  IMPLICIT NONE
  integer :: ierr
  write(*,*) "movetype_exit:"
  DEALLOCATE(movetype_name,stat=ierr)
  DEALLOCATE(movetype_prob,stat=ierr)
  DEALLOCATE(movetype_ntry,stat=ierr)
  DEALLOCATE(movetype_nsuccess,stat=ierr)
  DEALLOCATE(movetype_i_success,stat=ierr)
  DEALLOCATE(movetype_i_try,stat=ierr)
  if (ensemble_exch) then
    DEALLOCATE(exch_tcomp,stat=ierr)
  endif
END SUBROUTINE movetype_exit

subroutine movetype_select(guess,select_type)
    use movetype, only: nmovetypes, movetype_prob
    implicit none
    logical :: check_movetype
    double precision :: guess, tprob
    integer :: itype, select_type
    check_movetype = .false.
    tprob = 0.0D0
    do itype = 1, nmovetypes
      tprob = tprob + movetype_prob(itype)
      if (guess < tprob) then
        check_movetype = .true.
        exit
      endif
    enddo
    if (.not. check_movetype) THEN
      write(*,*) " wrong check_movetype"
      STOP
    endif
    select_type = itype
    return
end subroutine movetype_select

subroutine movetype_run(itype)
  use movetype, only: movetype_name
  use pos, only: box
  use ints, only: nptot
  use coupling_pres, only: dvx, delta, expd
  implicit none
  integer :: itype, imol
  double precision :: rand
  external rand
  if (movetype_name(itype) == 'trans') then
    imol = INT( NPTOT*RAND() ) + 1 ! pick random particle, (0 <= RAND < 1)
    call try_conf(imol,itype)
  else if (movetype_name(itype) == 'press') then
    delta = (2.0D0*RAND()-1.0D0)*dvx   ! pick volume change direction or total volume
    expd  = (box(3)+delta)/box(3)      ! get expansion ratio on z-axis (use this factor for other axes)
    call try_pres(itype)
  else if (movetype_name(itype) == 'exch') then
    imol = INT( NPTOT*RAND() ) + 1 ! pick random particle, (0 <= RAND < 1)
    call try_exch(imol,itype)
  else
    write(*,*) " wrong argument itype"
    stop
  endif
  return
end subroutine movetype_run


text:
module mod_monolis_precond_diag
  use mod_monolis_prm
  use mod_monolis_com
  use mod_monolis_mat
  use mod_monolis_precond_diag_33
  use mod_monolis_precond_diag_nn

  implicit none

contains

  subroutine monolis_precond_diag_setup(monoPRM, monoCOM, monoMAT)
    implicit none
    type(monolis_prm) :: monoPRM
    type(monolis_com) :: monoCOM
    type(monolis_mat) :: monoMAT

    if(monoMAT%NDOF == 3)then
      call monolis_precond_diag_33_setup(monoMAT)
    else
      call monolis_precond_diag_nn_setup(monoMAT)
    endif
  end subroutine monolis_precond_diag_setup

  subroutine monolis_precond_diag_apply(monoPRM, monoCOM, monoMAT, X, Y)
    implicit none
    type(monolis_prm) :: monoPRM
    type(monolis_com) :: monoCOM
    type(monolis_mat) :: monoMAT
    real(kind=kdouble) :: X(:), Y(:)

    if(monoMAT%NDOF == 3)then
      call monolis_precond_diag_33_apply(monoMAT, X, Y)
    else
      call monolis_precond_diag_nn_apply(monoMAT, X, Y)
    endif
  end subroutine monolis_precond_diag_apply

  subroutine monolis_precond_diag_clear(monoPRM, monoCOM, monoMAT)
    implicit none
    type(monolis_prm) :: monoPRM
    type(monolis_com) :: monoCOM
    type(monolis_mat) :: monoMAT

    if(monoMAT%NDOF == 3)then
      call monolis_precond_diag_33_clear(monoMAT)
    else
      call monolis_precond_diag_nn_clear(monoMAT)
    endif
  end subroutine monolis_precond_diag_clear
end module mod_monolis_precond_diag

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis


      SUBROUTINE DSINFT(Y,N)
! Double-precision version of Numerical Recipes fast sine transform.
      implicit none
      integer n, j
      real*8 y(n)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA,y1,y2,sum

      THETA=3.141592653589793D0/DBLE(N)
      WR=1.0D0
      WI=0.0D0
      WPR=-2.0D0*DSIN(0.5D0*THETA)**2
      WPI=DSIN(THETA)
      Y(1)=0.0d0
      DO 11 J=1,N/2
        WTEMP=WR
        WR=WR*WPR-WI*WPI+WR
        WI=WI*WPR+WTEMP*WPI+WI
        Y1=WI*(Y(J+1)+Y(N-J+1))
        Y2=0.5d0*(Y(J+1)-Y(N-J+1))
        Y(J+1)=Y1+Y2
        Y(N-J+1)=Y1-Y2
11    CONTINUE
      CALL DREALFT(Y,N,+1)
      SUM=0.0d0
      Y(1)=0.5d0*Y(1)
      Y(2)=0.0d0
      DO 12 J=1,N-1,2
        SUM=SUM+Y(J)
        Y(J)=Y(J+1)
        Y(J+1)=SUM
12    CONTINUE
      RETURN
      END


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE CMDELB (VERB, INLINE, IFLD, INTYP, CFIELD, IFIELD,
     &                   IDELB, IELBST, NEWELB, *)
C=======================================================================

C   --*** CMDELB *** (MESH) Process active element block commands
C   --   Written by Amy Gilkey - revised 05/18/88
C   --
C   --Parameters:
C   --   VERB    - I/O - the verbs for the SHOW command
C   --   INLINE  - I/O - the parsed input line for the log file
C   --   IFLD    - I/O - the free-field reader index and fields
C   --   INTYP   - I/O - the free-field reader index and fields
C   --   CFIELD  - I/O - the free-field reader index and fields
C   --   IFIELD  - I/O - the free-field reader index and fields
C   --   IDELB   - IN  - the element block IDs
C   --   IELBST  - I/O - the element block status:
C   --                   -1 = OFF, 0 = ON, but not selected, 1 = selected
C   --   NEWELB  - I/O - the new element blocks flag, set 0 in ROUTER only
C   --                   0 = no new element blocks
C   --                   1 = new selected element blocks
C   --                   2 = new displayed element blocks
C   --                       (implies new selected blocks)
C   --
C   --Common Variables:
C   --   Uses IS3DIM of /D3NUMS/
C   --   Uses NELBLK of /DBNUMS/

      include 'params.blk'
      include 'dbnums.blk'
      include 'd3nums.blk'

      CHARACTER*(*) VERB
      CHARACTER*(*) INLINE
      INTEGER     INTYP(*)
      CHARACTER*(*) CFIELD(*)
      INTEGER     IFIELD(*)
      INTEGER IDELB(NELBLK)
      INTEGER IELBST(NELBLK)
      INTEGER NEWELB

      CHARACTER*(MXSTLN) WORD
      CHARACTER OPTION
      LOGICAL FFEXST, FFNUMB, FFMATC

      IF (VERB .EQ. 'VISIBLE') THEN
         CALL FFADDC (VERB, INLINE)
         IF (.NOT. IS3DIM) THEN
            CALL PRTERR ('CMDERR', 'Command allowed in 3D only')
            GOTO 190
         END IF

         NEWELB = 2

         IF (.NOT. FFEXST (IFLD, INTYP)) THEN

C         --Select all element blocks if no parameters
            DO 100 IELB = 1, NELBLK
               IF (IELBST(IELB) .LT. 0) IELBST(IELB) = 1
  100       CONTINUE

         ELSE IF (FFMATC (IFLD, INTYP, CFIELD, 'OFF', 3)) THEN

C         --Turn off all element blocks if OFF
            CALL FFADDC ('OFF', INLINE)
            DO 110 IELB = 1, NELBLK
               IF (IELBST(IELB) .GE. 0) IELBST(IELB) = -1
  110       CONTINUE

         ELSE

C         --Strip off ADD or DELETE option
            IF (FFMATC (IFLD, INTYP, CFIELD, 'ADD', 3)) THEN
               CALL FFADDC ('ADD', INLINE)
               OPTION = '+'
            ELSE IF (FFMATC (IFLD, INTYP, CFIELD, 'DELETE', 3)) THEN
               CALL FFADDC ('DELETE', INLINE)
               OPTION = '-'
            ELSE
               IF (.NOT. FFNUMB (IFLD, INTYP)) THEN
                  CALL PRTERR ('CMDERR',
     &               'Expected "OFF" or "ADD" or "DELETE"'
     &               // ' or element block ID')
                  GOTO 190
               END IF
C            --De-select all element blocks so only listed blocks are selected
               CALL INIINT (NELBLK, -1, IELBST)
               OPTION = '+'
            END IF

  120       CONTINUE
            IF (FFEXST (IFLD, INTYP)) THEN
               CALL FFINTG (IFLD, INTYP, IFIELD,
     &            'element block ID', 0, ID, *130)
               IELB = LOCINT (ID, NELBLK, IDELB)
               IF (IELB .GT. 0) THEN
                  CALL FFADDI (ID, INLINE)
                  IF (OPTION .EQ. '+') THEN
                     IF (IELBST(IELB) .LT. 0) IELBST(IELB) = 1
                  ELSE
                     IELBST(IELB) = -1
                  END IF
               ELSE
                  CALL INTSTR (1, 0, ID, WORD, LSTR)
                  CALL PRTERR ('CMDWARN', 'Element block ID '
     &               // WORD(:LSTR) // ' does not exist, ignored')
               END IF
  130          CONTINUE
               GOTO 120
            END IF

            CALL CNTELB (IELBST, NELBLK, NUMON, NUMSEL)
            IF (NUMON .EQ. 0) CALL FFADDC ('OFF', INLINE)
         END IF

      ELSE IF ((VERB .EQ. 'BLOCKS') .OR. (VERB .EQ. 'MATERIAL')) THEN
         CALL FFADDC (VERB, INLINE)
         NEWELB = MAX (NEWELB, 1)

         IF (.NOT. FFEXST (IFLD, INTYP)) THEN

C         --Select all element blocks if no parameters
            DO 140 IELB = 1, NELBLK
               IF (IELBST(IELB) .GE. 0) IELBST(IELB) = 1
  140       CONTINUE

         ELSE IF (FFMATC (IFLD, INTYP, CFIELD, 'OFF', 3)) THEN

C         --De-select all element blocks if OFF
            CALL FFADDC ('OFF', INLINE)
            DO 150 IELB = 1, NELBLK
               IF (IELBST(IELB) .GE. 0) IELBST(IELB) = 0
  150       CONTINUE

         ELSE

C         --Strip off ADD or DELETE option
            IF (FFMATC (IFLD, INTYP, CFIELD, 'ADD', 3)) THEN
               CALL FFADDC ('ADD', INLINE)
               OPTION = '+'
            ELSE IF (FFMATC (IFLD, INTYP, CFIELD, 'DELETE', 3)) THEN
               CALL FFADDC ('DELETE', INLINE)
               OPTION = '-'
            ELSE
               IF (.NOT. FFNUMB (IFLD, INTYP)) THEN
                  CALL PRTERR ('CMDERR',
     &               'Expected "OFF" or "ADD" or "DELETE"'
     &               // ' or element block ID')
                  GOTO 190
               END IF
C            --De-select all element blocks so only listed blocks are selected
               DO 160 IELB = 1, NELBLK
                  IF (IELBST(IELB) .GE. 0) IELBST(IELB) = 0
  160          CONTINUE
               OPTION = '+'
            END IF

  170       CONTINUE
            IF (FFEXST (IFLD, INTYP)) THEN
               CALL FFINTG (IFLD, INTYP, IFIELD,
     &            'element block ID', 0, ID, *180)
               IELB = LOCINT (ID, NELBLK, IDELB)
               IF (IELB .GT. 0) THEN
                  IF (IELBST(IELB) .GE. 0) THEN
                     CALL FFADDI (ID, INLINE)
                     IF (OPTION .EQ. '+') THEN
                        IELBST(IELB) = 1
                     ELSE
                        IELBST(IELB) = 0
                     END IF
                  ELSE
                     CALL INTSTR (1, 0, ID, WORD, LSTR)
                     CALL PRTERR ('CMDWARN', 'Element block '
     &                  // WORD(:LSTR) // ' is not displayed, ignored')
                  END IF
               ELSE
                  CALL INTSTR (1, 0, ID, WORD, LSTR)
                  CALL PRTERR ('CMDWARN', 'Element block ID '
     &               // WORD(:LSTR) // ' does not exist, ignored')
               END IF
  180          CONTINUE
               GOTO 170
            END IF

            CALL CNTELB (IELBST, NELBLK, NUMON, NUMSEL)
            IF (NUMSEL .EQ. 0) CALL FFADDC ('OFF', INLINE)
         END IF
      END IF

      RETURN

  190 CONTINUE
      RETURN 1
      END

ext: f

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
c EQUISUBMU_contmag evalua la presion en equilibrio hidrostatico
c teniendo en cuenta la presion magnetica
c_______________________________________________________________
c equisubmu rutina que evalua la presion en equilibrio hidrostatico
c teniendo en cuenta el numero de electrones para el calculo de mu
c Luis Bellot y Basilio Ruiz 27/7/95
c Cambio del metodo de integracion Basilio Ruiz 3/9/96
c _______________________________________________________________
	subroutine equisubmu_contmag(ntau,tau1,t,pe,pg0,pg,z,ro,b,gamma)

	implicit real*4 (a-h,o-z)

	include 'PARAMETER'
	parameter (nex=28,cgases=83145100.)
	real*4 tau(kt),t(kt),pe(kt),pg(*),kac,d2,x(kt),kappa(kt),taue(kt)
        real*4 z(kt),z1(kt),ro(*),y(kt),b(kt),gamma(*),ptot(kt)
        real*4 tau1(kt)
	real*4 mu
	integer*4 nmaxitera

	real*4 wgt,abu,ei1,ei2,pp(10),tsi,psi,d1(10)
        common/constantes/g,avog	!gravedad,n. avogadro/pmu
	common/mu/cth                   !esto esta deshabilitado desde sir (entra 1)
        common/precisoitera/precitera
        common/anguloheliocent/mu
        common/nmaxitera/nmaxitera

        precitera=1.e-5
        nmaxitera=250

	g=mu*2.7414e+4		!gravedad cm/s^2 en fotosfera solar
       	avog=6.023e23
        epsilon=.95              !cota  epsilon=max(Pmag/Ptot)
        pi=3.1415916

        do i=1,ntau
           tau(i)=tau1(i)+alog10(abs(cth))  !this is the optical depth in
        end do                        !the VERTICAL direction for
	                              !hydrostatic equilib. computations

        tau0=1.
        imin=1!esto esta deshabilitado desde sir (entra 1)
        do i=2,ntau
           if(abs(tau(i)).lt.tau0)then
              imin=i         !zero of z is at log tau (vertical) =0
              tau0=tau(i)
           end if
        end do

        do i=1,10
           d1(i)=0
        end do
c calculamos el peso molecular medio pmu
	pmusum=0.0
        asum=0.0
	do i=1,92
  	   ii=i
	   call neldatb(ii,0.,wgt,abu,ei1,ei2)
	   pmusum=pmusum+wgt*abu
           asum=asum+abu
	end do

	do i=1,ntau
           taue(i)=10.**tau(i)
        end do

	do i=1,ntau-1
           x(i+1)=g*(tau(i)-tau(i+1))*2.3025851
        end do

c inicializamos (atencion gamma debe entrar en radianes)
	tsi=t(ntau)
        psi=pe(ntau)
        pmagsi=b(ntau)*b(ntau)/8./pi*sin(gamma(ntau))*sin(gamma(ntau))
        if(pmagsi.gt.epsilon*pg0)pmagsi=epsilon*pg0

        pg00=pg0-pmagsi !nueva presion gaseosa
c        print*,'pg0 pg00 pmagsi epsilon',pg0,pg00,pmagsi,epsilon,sin(gamma(ntau))

        call pefrompg11(tsi,pg00,psi)
        pe(ntau)=psi
        call gasc(tsi,psi,pg(ntau),pp)
        pg(ntau)=pg00
        ptot(ntau)=pg0
	call kappach(5.e-5,tsi,psi,pp,d1,d1,kac,d2,d2)
        pesomedio=pmusum/(asum+pp(8)) !peso molec. medio
        kappa(ntau)=kac*avog/pmusum
        ro(ntau)=pesomedio*pg(ntau)/tsi/cgases
        y(ntau)=taue(ntau)/kappa(ntau)/ro(ntau)

c integramos
        do i=ntau-1,1,-1
           ptotpr=ptot(i+1)+x(i+1)*(taue(i+1)/kappa(i+1))
	   tsi=t(i)

           pmagsi=b(i)*b(i)/8./pi*sin(gamma(i))*sin(gamma(i))
           if(pmagsi.gt.epsilon*ptotpr)pmagsi=epsilon*ptotpr
           pgpr=ptotpr-pmagsi
           call pefrompg11(tsi,pgpr,psi)
           pe(i)=psi

           call gasc(tsi,psi,pgpr,pp)
	   call kappach(5.e-5,tsi,psi,pp,d1,d1,kac,d2,d2)
           kappa(i)=kac*avog/pmusum

           ptot(i)=ptot(i+1)+x(i+1)*(taue(i+1)/kappa(i+1)+taue(i)/kappa(i))/2.

           nit=1
           dif=2.*abs((ptot(i)-ptotpr))/(ptot(i)+ptotpr)

           do while (nit.lt.30.and.dif.gt.7.e-6)
              ptotpr=ptot(i)
              nit=nit+1
              pmagsi=b(i)*b(i)/8./pi*sin(gamma(i))*sin(gamma(i))
              if(pmagsi.gt.epsilon*ptotpr)pmagsi=epsilon*ptotpr
              pgpr=ptotpr-pmagsi

              call pefrompg11(tsi,pgpr,psi)
              pe(i)=psi
              call gasc(tsi,psi,pgpr,pp)
	      call kappach(5.e-5,tsi,psi,pp,d1,d1,kac,d2,d2)
              kappa(i)=kac*avog/pmusum
            ptot(i)=ptot(i+1)+x(i+1)*(taue(i+1)/kappa(i+1)+taue(i)/kappa(i) )/2.
              dif=2.*abs((ptot(i)-ptotpr))/(ptot(i)+ptotpr)
            end do

            if(dif.gt.0.1) print*,'WARNING:
     & Hydrostatic equilibrium results in inaccurate electron pressures '
            pg(i)=pgpr
            call pefrompg11(tsi,pg(i),psi)
            pe(i)=psi
            ro(i)=pesomedio*pg(i)/tsi/cgases
            y(i)=taue(i)/kappa(i)/ro(i)
	end do
        z1(1)=0.
        do i=2,ntau
           z1(i)=z1(i-1)+x(i)/g*(y(i-1)+y(i))/2.
        end do
	z00=z1(imin)

        do i=1,ntau
           z(i)=(z1(i)-z00)*1.e-5
        end do

c        do i=1,ntau   !we output geometrical heights ALONG the LOS
c          z(i)=z(i)/abs(cth)
c        enddo

       return
       end

max_stars_repo_name: cdiazbas/SIRcode

max_issues_repo_name: cdiazbas/SIRcode

max_forks_repo_name: cdiazbas/SIRcode


text:
	SUBROUTINE IP_GTPT  ( pvar, iret )
C************************************************************************
C* IP_GTPT								*
C*									*
C* This subroutine gets the points from the button presses and sets	*
C* the requested variable.						*
C*									*
C* IP_GTPT  ( PVAR, IRET )						*
C*									*
C* Input parameters:							*
C*	PVAR		CHAR*		Input parameter name		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*				   	  0 = normal return 		*
C*					-11 = Invalid parameter name	*
C*					-41 = Invalid cursor points	*
C**									*
C* Log:									*
C* S. Jacobs/EAI	 6/93						*
C* S. Jacobs/EAI	 9/93		Added ITYP to call for GGTPNT	*
C* S. Jacobs/NMC	 6/94		Changed values for parms to	*
C*					one decimal place		*
C* K. Tyle/GSC		 7/96		Renamed from NT_GTPT		*
C* K. Tyle/GSC		 7/96		Added check for mode		*
C* S. Maxwell/GSC        1/97           Added ST_LSTF, ST_LSTC, and     *
C*                                      GG_ZARE                         *
C* S. Jacobs/NCEP	 1/97		Fixed check for SAT projection;	*
C*					Removed unused variables	*
C* D.W.Plummer/NCEP	 6/97	Increased string lengths from 72 to 128	*
C* S. Jacobs/NCEP	 6/98		Removed NP from call to GGTPNT	*
C* T. Lee/GSC		 1/01		Added extended display area	*
C* M. Li/SAIC		12/07		Add CGAREA and CGAREAX		*
C************************************************************************
	INCLUDE		'ipcmn.cmn'
C*
	CHARACTER*(*)	pvar
C*
	REAL		x(2), y(2), rarr(2)
	CHARACTER	value*128, wname*128, garea*128
	CHARACTER	proj*128, str(2)*128
C-----------------------------------------------------------------------
	iret = 0
	CALL ST_LCUC ( pvar, pvar, ier )
C
C*	Get the integer from the input.
C
	IF  ( pvar .eq. 'GPOINT' )  THEN
	    ityp = 1
	ELSE IF  ( pvar .eq. 'CXSTNS' )  THEN
	    ityp = 2
	ELSE IF  ( pvar .eq. 'GAREA' )  THEN
	    ityp = 3
	ELSE IF  ( pvar .eq. 'GAREAX' )  THEN
	    ityp = 13
	    pvar = 'GAREA'
	ELSE IF  ( pvar .eq. 'CGAREA' )  THEN
            ityp = 4
            pvar = 'GAREA'
	ELSE IF  ( pvar .eq. 'CGAREAX' )  THEN
            ityp = 14
            pvar = 'GAREA'
	ELSE
	    iret = -11
	    CALL ER_WMSG ( 'IP', iret, ' ', ier )
	    RETURN
	END IF
C
C*	Select and pop the current window.
C
	wname = ' '
	CALL GSLWIN ( wname, ier )
	CALL GEPLOT ( ier )
C
C*	Get the cursor point(s).
C
	CALL GGTPNT ( 'M', ityp, x, y, iret )
	IF  ( iret .ne. 0 )  RETURN
C
C*	Set the value of the parameter.
C
	IF  ( pvar .eq. 'GPOINT' )  THEN
       	    rarr(1) = x(1)
            rarr(2) = y(1)
            CALL ST_LSTF ( rarr, 2, ';', 1, value, iret )
	  ELSE IF  ( pvar .eq. 'CXSTNS' )  THEN
            rarr(1) = x(1)
            rarr(2) = y(1)
            CALL ST_LSTF ( rarr, 2, ';', 1, str(1), iret )
            rarr(1) = x(2)
            rarr(2) = y(2)
            CALL ST_LSTF ( rarr, 2, ';', 1, str(2), iret )
            CALL ST_LSTC ( str, 2, '>', value, iret )
  	  ELSE IF  ( pvar .eq. 'GAREA' )  THEN
C
C*	    Check for map mode.
C
	    CALL GQMODE ( mode, iret )
	    IF ( mode .eq. 1 ) THEN
C
C*		Get new garea and proj.
C
		CALL GG_ZARE ( x, y, garea, proj, iret )
		value = garea
C
C*	   	Check value of PROJ.
C
	    	DO i = 1, ncparm
		    IF ( cparmn (i) .eq. 'PROJ' ) ilist = i
	    	END DO
C
	    	CALL ST_LCUC ( cparmv (ilist), cparmv (ilist), ier)
C
C*		Redefine value of PROJ, if it is not currently
C*		SAT or RAD.
C
	    	IF  ( ( cparmv (ilist) .ne. 'SAT' ) .and.
     +		      ( cparmv (ilist) .ne. 'RAD') ) THEN
		    cparmv (ilist) = proj
	    	END IF
	    END IF
	END IF
C
C*	Check for parameter name in full list.
C
	DO  i = 1, ncparm
	    IF  ( pvar .eq. cparmn (i) )  ifound = i
	END DO
	cparmv ( ifound ) = value
C
C*	Flush and pop the current window.
C
	CALL GEPLOT ( ier )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      subroutine clawpack46_bc2_default(maxmx,maxmy,meqn,mbc,
     &      mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt,mthbc)
      implicit none

      integer maxmx, maxmy, meqn, mbc, mx, my, maux, mthbc(4)
      double precision xlower, ylower, dx, dy, t, dt

      double precision q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, meqn)
      double precision aux(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, maux)

      integer m, i, j, ibc, jbc

c
c-------------------------------------------------------
c     # left boundary:
c-------------------------------------------------------
      go to (100,110,120,130) mthbc(1)+1
c     this is how we skip over this side... if (mthbc(1)+1
c     is not 1,2,3 or 4, then the goto above falls through to here...
      goto 199

  100 continue
c     # user-specified boundary conditions go here in place of error output
      write(6,*) '*** ERROR *** mthbc(1)=0 and no BCs specified in bc2'
      stop
      go to 199
c
  110 continue
c     # zero-order extrapolation:
      do m=1,meqn
         do ibc=1,mbc
            do j = 1-mbc,my+mbc
               q(1-ibc,j,m) = q(1,j,m)
            end do
         end do
      end do
      go to 199

  120 continue
c     # periodic:
      write(6,*) 'clawpack4.6_bc2_default : ',
     &            'Periodic BCs should not be set in bc2'
      stop
      go to 199

  130 continue
c     # solid wall (assumes 2'nd component is velocity or momentum in x):
      do m=1,meqn
         do ibc=1,mbc
            do j = 1-mbc, my+mbc
               q(1-ibc,j,m) = q(ibc,j,m)
            end do
         end do
      end do
c     # negate the normal velocity:
      do ibc=1,mbc
         do j = 1-mbc, my+mbc
            q(1-ibc,j,2) = -q(ibc,j,2)
         end do
      end do
      go to 199

  199 continue
c
c-------------------------------------------------------
c     # right boundary:
c-------------------------------------------------------
      go to (200,210,220,230) mthbc(2)+1
      goto 299
c
  200 continue
c     # user-specified boundary conditions go here in place of error output
      write(6,*) '*** ERROR *** mthbc(2)=0 and no BCs specified in bc2'
      stop
      go to 299

  210 continue
c     # zero-order extrapolation:
      do m=1,meqn
         do ibc=1,mbc
            do j = 1-mbc,my+mbc
               q(mx+ibc,j,m) = q(mx,j,m)
            end do
         end do
      end do
      go to 299

  220 continue
c     # periodic:
      write(6,*) 'clawpack4.6_bc2_default : ',
     &            'Periodic BCs should not be set in bc2'
      stop
      go to 299

  230 continue
c     # solid wall (assumes 2'nd component is velocity or momentum in x):
      do m=1,meqn
         do ibc=1,mbc
            do j = 1-mbc, my+mbc
               q(mx+ibc,j,m) = q(mx+1-ibc,j,m)
            end do
         end do
      end do
c     # negate the normal velocity:
      do ibc=1,mbc
         do j = 1-mbc, my+mbc
            q(mx+ibc,j,2) = -q(mx+1-ibc,j,2)
         end do
      end do
      go to 299

  299 continue
c
c-------------------------------------------------------
c     # bottom boundary:
c-------------------------------------------------------
      go to (300,310,320,330) mthbc(3)+1
      goto 399
c
  300 continue
c     # user-specified boundary conditions go here in place of error output
      write(6,*) '*** ERROR *** mthbc(3)=0 and no BCs specified in bc2'
      stop
      go to 399
c
  310 continue
c     # zero-order extrapolation:
      do m=1,meqn
         do jbc=1,mbc
            do i = 1-mbc,mx+mbc
               q(i,1-jbc,m) = q(i,1,m)
            end do
         end do
      end do
      go to 399

  320 continue
c     # periodic:
      write(6,*) 'clawpack4.6_bc2_default : ',
     &            'Periodic BCs should not be set in bc2'
      stop
      go to 399

  330 continue
c     # solid wall (assumes 3'rd component is velocity or momentum in y):
      do m=1,meqn
         do jbc=1,mbc
            do i = 1-mbc,mx+mbc
               q(i,1-jbc,m) = q(i,jbc,m)
            end do
         end do
      end do
c     # negate the normal velocity:
      do jbc=1,mbc
         do i = 1-mbc, mx+mbc
            q(i,1-jbc,3) = -q(i,jbc,3)
         end do
      end do
      go to 399

  399 continue
c
c-------------------------------------------------------
c     # top boundary:
c-------------------------------------------------------
      go to (400,410,420,430) mthbc(4)+1
      goto 499
c
  400 continue
c     # user-specified boundary conditions go here in place of error output
      write(6,*) '*** ERROR *** mthbc(4)=0 and no BCs specified in bc2'
      stop
      go to 499

  410 continue
c     # zero-order extrapolation:
      do m=1,meqn
         do jbc=1,mbc
            do i = 1-mbc,mx+mbc
               q(i,my+jbc,m) = q(i,my,m)
            end do
         end do
      end do
      go to 499

  420 continue
c     # periodic:
      write(6,*) 'clawpack4.6_bc2_default : ',
     &            'Periodic BCs should not be set in bc2'
      stop
      go to 499

  430 continue
c     # solid wall (assumes 3'rd component is velocity or momentum in y):
      do m=1,meqn
         do jbc=1,mbc
            do i = 1-mbc, mx+mbc
               q(i,my+jbc,m) = q(i,my+1-jbc,m)
            end do
         end do
      end do
c     # negate the normal velocity:
      do jbc=1,mbc
         do i = 1-mbc, mx+mbc
            q(i,my+jbc,3) = -q(i,my+1-jbc,3)
         end do
      end do
      go to 499

  499 continue

      return
      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


!> \file mo_linear_algebra.f90

!> \brief Wrapper for LAPACK's linear algebra routines.

!> \details This modules provides mostly wrappers for LAPACK's F77 linear algebra routines.
!>          It adds a few convenience functions such as diag.

!> \authors Matthias Cuntz, Sebastian Mueller
!> \date May 2014
MODULE mo_linear_algebra

  ! Wrapper for LAPACK's F77 linear algebra routines.

  ! Written  Matthias Cuntz,    May 2014
  ! Modified Matthias Cuntz,    May 2016 - calc single precision via double precision
  ! Modified Sebastian Mueller, Oct 2016 - solver for banded coefficent-matrices and some involved algorithms
  ! Modified Matthias Cuntz,    Mar 2020 - allocate out of sp routines calling dp routines only if not allocated
  !                                      - allocate out for Python.

  ! License
  ! -------
  ! This file is part of the JAMS Fortran package, distributed under the MIT License.
  !
  ! Copyright (c) 2014-2020 Matthias Cuntz, Sebastian Mueller - mc (at) macu (dot) de
  !
  ! 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.

  USE mo_kind, ONLY: i4, i8, sp, dp

  IMPLICIT NONE

  PUBLIC :: banded                      ! banded form of a given matrix
  PUBLIC :: diag                        ! diagonal of matrix
  PUBLIC :: inverse                     ! inverse of matrix
  PUBLIC :: min_diag                    ! minor diagonal of matrix
  PUBLIC :: solve_linear_equations      ! solve linear system of equations with LU decomposition
  PUBLIC :: solve_linear_equations_svd  ! solve linear system of equations with SVD
  PUBLIC :: solve_linear_equations_band ! solve linear system of equations for banded matrix

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         banded
  !
  !     PURPOSE
  !         Banded form of squared matrix for the banded linear equations solver
  !
  !>        \brief Banded form of squared matrix.
  !
  !>        \details Converts a given squared matrix into the banded form needed for the Lapack-solver for banded matrices.
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp) :: matrix(:,:)"    Squared 2D-array
  !>        \param[in] "integer(i4) :: l          "    number of lower minor-diagonals
  !>        \param[in] "integer(i4) :: u          "    number of upper minor-diagonals
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !         None
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp) :: banded(:) — Banded form of input matrix
  !
  !     RESTRICTIONS
  !         None
  !
  !     EXAMPLE
  !         band = banded(matrix,l,u)
  !         -> see also example in test directory

  !     LITERATURE
  !         None

  !     HISTORY
  !>        \author Sebastian Mueller
  !>        \date October 2016
  INTERFACE banded
     MODULE PROCEDURE banded_dp, banded_sp
  END INTERFACE banded

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         diag
  !
  !     PURPOSE
  !         Returns the diagonal of a square 2D-array.
  !
  !>        \brief Diagonal elements of a squared matrix.
  !
  !>        \details Returns the diagonal of a square 2D-array.
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp)/integer(i4/i8)/logical :: matrix(:,:)"    Squared 2D-array
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !         None
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp)/integer(i4/i8)/logical :: diag(:) — Diagonal elements of matrix
  !
  !     RESTRICTIONS
  !         None
  !
  !     EXAMPLE
  !         m = diag(matrix)
  !         -> see also example in test directory
  !
  !     LITERATURE
  !         None
  !
  !     HISTORY
  !>        \author Matthias Cuntz
  !>        \date May 2014
  INTERFACE diag
     MODULE PROCEDURE diag_sp, diag_dp, diag_i4, diag_i8, diag_lgt
  END INTERFACE diag

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         inverse
  !
  !     PURPOSE
  !         Inverse of squared matrix
  !
  !>        \brief Inverse of squared matrix.
  !
  !>        \details Inverts squared matrix using LU decomposition.
  !>
  !>                 Uses standard Lapack routines using LU decomposition: dgetrf and dgetri.
  !>                 Conditions columns before decomposition.
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp) :: matrix(:,:)"    Squared 2D-array
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !>       \param[in] "logical, optional :: condition" If true, condition matrix before inversion (default: true)
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp) :: inverse(:) — Inverse of input matrix
  !
  !     RESTRICTIONS
  !         None
  !
  !     EXAMPLE
  !         inv = inverse(matrix)
  !         -> see also example in test directory

  !     LITERATURE
  !         None

  !     HISTORY
  !>        \author Matthias Cuntz
  !>        \date May 2014
  INTERFACE inverse
     MODULE PROCEDURE inverse_dp, inverse_sp
  END INTERFACE inverse

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         min_diag
  !
  !     PURPOSE
  !         Returns the n-th minor-diagonal of a square 2D-array, where positive n indicate super-diagonals and negativ n indicate
  !         sub-diagonals.
  !
  !>        \brief Minor-diagonal elements of a squared matrix.
  !
  !>        \details Returns the n-th minor-diagonal of a square 2D-array.
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp)/integer(i4/i8)/logical :: matrix(:,:)"    Squared 2D-array
  !>        \param[in] "integer(i4)                        :: n          "    number of minor-diagonal
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !         None
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp)/integer(i4/i8)/logical :: min_diag(:) — n-th minor-diagonal elements of matrix
  !
  !     RESTRICTIONS
  !         None
  !
  !     EXAMPLE
  !         m = min_diag(matrix,n)
  !         -> see also example in test directory
  !
  !     LITERATURE
  !         None
  !
  !     HISTORY
  !>        \author Sebastian Mueller
  !>        \date October 2016
  INTERFACE min_diag
     MODULE PROCEDURE min_diag_sp, min_diag_dp, min_diag_i4, min_diag_i8, min_diag_lgt
  END INTERFACE min_diag

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         solve_linear_equations
  !
  !     PURPOSE
  !         Solve linear system of equations
  !
  !>        \brief Solve linear system of equations.
  !
  !>        \details Solve linear system of equations using LU decomposition
  !>                 \f[ A x = b \f]
  !>                 Returns \f$ x \f$.
  !>
  !>                 Uses standard Lapack routine using LU decomposition: dgesv.
  !>                 Conditions rows and columns before decomposition.
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp) :: lhs(:,:)"    Coefficients of left hand side \f$ A \f$
  !>        \param[in] "real(sp/dp) :: rhs(:)"      Right hand side \f$ b \f$
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !>       \param[in] "logical, optional :: condition" If true, condition lhs and rhs before decomposition (default: true)
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp) :: x(:) — Solution x to \f$ A x = b \f$
  !
  !     RESTRICTIONS
  !         Only one right-hand side.
  !
  !     EXAMPLE
  !         sol = solve_linear_equations(lhs, rhs)
  !         -> see also example in test directory

  !     LITERATURE
  !         None

  !     HISTORY
  !>        \author Matthias Cuntz
  !>        \date May 2014
  INTERFACE solve_linear_equations
     MODULE PROCEDURE solve_linear_equations_1_dp, solve_linear_equations_1_sp
  END INTERFACE solve_linear_equations

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         solve_linear_equations_svd
  !
  !     PURPOSE
  !         Solve linear system of equations
  !
  !>        \brief Solve linear system of equations.
  !
  !>        \details Solve linear system of equations using singular value decomposition
  !>                 \f[ A x = b \f]
  !>                 Returns \f$ x \f$.
  !>
  !>                 Uses standard Lapack routine to do SVD: dgesvd.
  !>                 then solving routine similar to Numerical Recipes: svdksb.
  !>                 It conditions rows and columns before decomposition.
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp) :: lhs(:,:)"    Coefficients of left hand side \f$ A \f$
  !>        \param[in] "real(sp/dp) :: rhs(:)"      Right hand side \f$ b \f$
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !>       \param[in] "logical, optional :: condition" If true, condition lhs and rhs before decomposition (default: true)
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp) :: x(:) — Solution x to \f$ A x = b \f$
  !
  !     RESTRICTIONS
  !         Only one right-hand side.
  !
  !     EXAMPLE
  !         sol = solve_linear_equations_svd(lhs, rhs)
  !         -> see also example in test directory

  !     LITERATURE
  !         Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 -
  !             The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes,
  !             Cambridge University Press, UK, 1996

  !     HISTORY
  !>        \author Matthias Cuntz
  !>        \date May 2014
  INTERFACE solve_linear_equations_svd
     MODULE PROCEDURE solve_linear_equations_svd_1_dp, solve_linear_equations_svd_1_sp
  END INTERFACE solve_linear_equations_svd

  ! ------------------------------------------------------------------
  !
  !     NAME
  !         solve_linear_equations_band
  !
  !     PURPOSE
  !         Solve linear system of equations for a banded matrix A=(aij).
  !         The band storage scheme is illustrated by the following example, when
  !         N = 6, l = 2, u = 1:
  !
  !             *   a12  a23  a34  a45  a56
  !            a11  a22  a33  a44  a55  a66
  !            a21  a32  a43  a54  a65   *
  !            a31  a42  a53  a64   *    *
  !
  !         Array elements marked * are not used by the routine. You can set them zero.
  !
  !>        \brief Solve linear system of equations.
  !
  !>        \details Solve linear system of equations for banded coefficent-matrix
  !>                 \f[ A x = b \f]
  !>                 Returns \f$ x \f$.
  !>
  !>                 Uses standard Lapack routine for banded matrix: dgbsv.
  !>                 Conditions columns before decomposition.
  !
  !     INTENT(IN)
  !>        \param[in] "real(sp/dp) :: lhsb(:,:)"    Coefficients of left hand side in banded form of \f$ A \f$
  !>        \param[in] "real(sp/dp) :: rhs(:)   "    Right hand side \f$ b \f$
  !>        \param[in] "integer(i4) :: l        "    number of lower minor-diagonals
  !>        \param[in] "integer(i4) :: u        "    number of upper minor-diagonals
  !
  !     INTENT(INOUT)
  !         None
  !
  !     INTENT(OUT)
  !         None
  !
  !     INTENT(IN), OPTIONAL
  !>       \param[in] "logical, optional :: condition" If true, condition lhsb before decomposition (default: true)
  !
  !     INTENT(INOUT), OPTIONAL
  !         None
  !
  !     INTENT(OUT), OPTIONAL
  !         None
  !
  !     RETURNS
  !>       \return     real(sp/dp) :: x(:) — Solution x to \f$ A x = b \f$
  !
  !     RESTRICTIONS
  !         Only one right-hand side.
  !         Left-hand side needs to be in banded form. You can use the banded-function to convert a given squared matrix.
  !
  !     EXAMPLE
  !         sol = solve_linear_equations_band(lhsb, rhs, l, u)
  !         sol = solve_linear_equations_band(banded(lhs), rhs, l, u)
  !         -> see also example in test directory

  !     LITERATURE
  !         None

  !     HISTORY
  !>        \author Sebastian Mueller
  !>        \date October 2016
  !         Modified, Matthias Cuntz, May 2018 - allocate output for solve_linear_equations_band_1_sp
  INTERFACE solve_linear_equations_band
     MODULE PROCEDURE solve_linear_equations_band_1_dp, solve_linear_equations_band_1_sp
  END INTERFACE solve_linear_equations_band

  PRIVATE

  INTERFACE svdksb
     MODULE PROCEDURE svdksb_dp, svdksb_sp
  END INTERFACE svdksb

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

CONTAINS

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

  FUNCTION banded_dp(matrix, l, u)

    IMPLICIT NONE

    REAL(dp), DIMENSION(:,:), INTENT(IN)   :: matrix
    INTEGER(i4),              INTENT(IN)   :: l, u
#ifndef __PYTHON__
    REAL(dp), DIMENSION(:,:), allocatable  :: banded_dp
#else
    REAL(dp), DIMENSION(l+u+1,size(matrix,1)) :: banded_dp
#endif

    INTEGER(i4) :: i

    if (size(matrix,1) /= size(matrix,2)) stop 'banded_dp: array must be squared matrix.'
    if (l              >= size(matrix,1)) stop 'banded_dp: l is to big. You need to choose minor diagonals within the matrix.'
    if (u              >= size(matrix,1)) stop 'banded_dp: u is to big. You need to choose minor diagonals within the matrix.'
    if (l              <  0_i4          ) stop 'banded_dp: l needs to be non-negativ.'
    if (u              <  0_i4          ) stop 'banded_dp: u needs to be non-negativ.'
#ifndef __PYTHON__
    if (.not. allocated(banded_dp)) allocate(banded_dp(l+u+1,size(matrix,1)))
#endif

    banded_dp = 0.0_dp

    do i=0_i4,u
        banded_dp(u+1-i,i+1:size(matrix,1)) = min_diag_dp(matrix, i)
    end do

    do i=1_i4,l
        banded_dp(u+1+i,1:size(matrix,1)-i) = min_diag_dp(matrix,-i)
    end do

  END FUNCTION banded_dp

  FUNCTION banded_sp(matrix, l, u)

    IMPLICIT NONE

    REAL(sp), DIMENSION(:,:), INTENT(IN)   :: matrix
    INTEGER(i4),              INTENT(IN)   :: l, u
#ifndef __PYTHON__
    REAL(sp), DIMENSION(:,:), allocatable  :: banded_sp
#else
    REAL(sp), DIMENSION(l+u+1,size(matrix,1)) :: banded_sp
#endif

    INTEGER(i4) :: i

    if (size(matrix,1) /= size(matrix,2)) stop 'banded_sp: array must be squared matrix.'
    if (l              >= size(matrix,1)) stop 'banded_sp: l is to big. You need to choose minor diagonals within the matrix.'
    if (u              >= size(matrix,1)) stop 'banded_sp: u is to big. You need to choose minor diagonals within the matrix.'
    if (l              <  0_i4          ) stop 'banded_sp: l needs to be non-negativ.'
    if (u              <  0_i4          ) stop 'banded_sp: u needs to be non-negativ.'
#ifndef __PYTHON__
    if (.not. allocated(banded_sp)) allocate(banded_sp(l+u+1,size(matrix,1)))
#endif

    banded_sp = 0.0_sp

    do i=0_i4,u
        banded_sp(u+1-i,i+1:size(matrix,1)) = min_diag_sp(matrix, i)
    end do

    do i=1_i4,l
        banded_sp(u+1+i,1:size(matrix,1)-i) = min_diag_sp(matrix,-i)
    end do

  END FUNCTION banded_sp

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

  FUNCTION diag_dp(matrix)

    IMPLICIT NONE

    REAL(dp), DIMENSION(:,:), INTENT(IN) :: matrix
#ifndef __PYTHON__
    REAL(dp), DIMENSION(:), allocatable  :: diag_dp
#else
    REAL(dp), DIMENSION(size(matrix,1)) :: diag_dp
#endif

    INTEGER(i4) :: i

    if (size(matrix,1) /= size(matrix,2)) stop 'diag_dp: array must be squared matrix.'
#ifndef __PYTHON__
    if (.not. allocated(diag_dp)) allocate(diag_dp(size(matrix,1)))
#endif

    forall(i=1:size(matrix,1)) diag_dp(i) = matrix(i,i)

  END FUNCTION diag_dp

  FUNCTION diag_sp(matrix)

    IMPLICIT NONE

    REAL(sp), DIMENSION(:,:), INTENT(IN) :: matrix
#ifndef __PYTHON__
    REAL(sp), DIMENSION(:), allocatable  :: diag_sp
#else
    REAL(sp), DIMENSION(size(matrix,1)) :: diag_sp
#endif

    INTEGER(i4) :: i

    if (size(matrix,1) /= size(matrix,2)) stop 'diag_sp: array must be squared matrix.'
#ifndef __PYTHON__
    if (.not. allocated(diag_sp)) allocate(diag_sp(size(matrix,1)))
#endif

    forall(i=1:size(matrix,1)) diag_sp(i) = matrix(i,i)

  END FUNCTION diag_sp

  FUNCTION diag_i4(matrix)

    IMPLICIT NONE

    INTEGER(i4), DIMENSION(:,:), INTENT(IN) :: matrix
#ifndef __PYTHON__
    INTEGER(i4), DIMENSION(:), allocatable  :: diag_i4
#else
    INTEGER(i4), DIMENSION(size(matrix,1)) :: diag_i4
#endif

    INTEGER(i4) :: i

    if (size(matrix,1) /= size(matrix,2)) stop 'diag_i4: array must be squared matrix.'
#ifndef __PYTHON__
    if (.not. allocated(diag_i4)) allocate(diag_i4(size(matrix,1)))
#endif

    forall(i=1:size(matrix,1)) diag_i4(i) = matrix(i,i)

  END FUNCTION diag_i4

  FUNCTION diag_i8(matrix)

    IMPLICIT NONE

    INTEGER(i8), DIMENSION(:,:), INTENT(IN) :: matrix
#ifndef __PYTHON__
    INTEGER(i8), DIMENSION(:), allocatable  :: diag_i8
#else
    INTEGER(i8), DIMENSION(size(matrix,1)) :: diag_i8
#endif

    INTEGER(i8) :: i

    if (size(matrix,1) /= size(matrix,2)) stop 'diag_i8: array must be squared matrix.'
#ifndef __PYTHON__
    if (.not. allocated(diag_i8)) allocate(diag_i8(size(matrix,1)))
#endif

    forall(i=1:size(matrix,1)) diag_i8(i) = matrix(i,i)

  END FUNCTION diag_i8

  FUNCTION diag_lgt(matrix)

    IMPLICIT NONE

    LOGICAL, DIMENSION(:,:), INTENT(IN) :: matrix
#ifndef __PYTHON__
    LOGICAL, DIMENSION(:), allocatable  :: diag_lgt
#else
    LOGICAL, DIMENSION(size(matrix,1)) :: diag_lgt
#endif

    INTEGER(i4) :: i

    if (size(matrix,1) /= size(matrix,2)) stop 'diag_lgt: array must be squared matrix.'
#ifndef __PYTHON__
    if (.not. allocated(diag_lgt)) allocate(diag_lgt(size(matrix,1)))
#endif

    forall(i=1:size(matrix,1)) diag_lgt(i) = matrix(i,i)

  END FUNCTION diag_lgt

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

  FUNCTION inverse_dp(matrix, condition)

    IMPLICIT NONE

    REAL(dp), DIMENSION(:,:), INTENT(IN)  :: matrix
    LOGICAL,  OPTIONAL,       INTENT(IN)  :: condition
#ifndef __PYTHON__
    REAL(dp), DIMENSION(:,:), allocatable :: inverse_dp
#else
    REAL(dp), DIMENSION(size(matrix,2),size(matrix,1)) :: inverse_dp
#endif

    INTEGER(i4) :: i, nn
    real(dp),    dimension(:),   allocatable :: scale_cols    ! scale matrix for better conditioning
    integer(i4), dimension(:),   allocatable :: ipiv          ! needed for dgetrf lapack routine
    integer(i4)                              :: info          !              "
    real(dp),    dimension(:),   allocatable :: work          ! needed for dgetri lapack routine
    integer(i4)                              :: lwork         !              "
    logical :: icondition

    external :: dgetrf ! Lapack routine to compute LU factorization of a general matrix
    external :: dgetri ! Lapack routine to compute inverse of matrix using the LU factorization computed by DGETRF

    nn = size(matrix,2)
    if (size(matrix,1) /= nn) stop 'inverse_dp: matrix must be square.'
#ifndef __PYTHON__
    if (.not. allocated(inverse_dp)) allocate(inverse_dp(nn,nn))
#endif
    inverse_dp = matrix

    if (present(condition)) then
       icondition = condition
    else
       icondition = .true.
    endif

    if (icondition) then
       ! Condition columns
       allocate(scale_cols(nn))
       scale_cols(:) = maxval(abs(inverse_dp(:,:)),1)
       where (scale_cols(:) < tiny(1.0_dp)) scale_cols(:) = 1.0_dp
       scale_cols(:) = 1.0_dp / scale_cols(:)
       forall(i=1:nn) inverse_dp(:,i) = inverse_dp(:,i) * scale_cols(i)
    endif

    ! LU factorisation of imatrix
    allocate(ipiv(nn))
    call dgetrf(nn, nn, inverse_dp, nn, ipiv, info)
    if (info /= 0) stop 'inverse_dp: LU factorisation did not work.'
    ! Inverse of LU factorisation of imatrix
    allocate(work(1))
    call dgetri(nn, inverse_dp, nn, ipiv, work, -1, info)
    lwork = int(work(1),i4)
    deallocate(work)
    allocate(work(lwork))
    call dgetri(nn, inverse_dp, nn, ipiv, work, lwork, info)
    if (info /= 0) stop 'hdmr_hessian: Inversion did not work.'

    if (icondition) then
       ! rescale result
       forall(i=1:nn) inverse_dp(i,:) = inverse_dp(i,:) * scale_cols(i)

       deallocate(scale_cols)
    endif

    deallocate(ipiv)

  END FUNCTION inverse_dp

  FUNCTION inverse_sp(matrix, condition)

    IMPLICIT NONE

    REAL(sp), DIMENSION(:,:), INTENT(IN)  :: matrix
    LOGICAL,  OPTIONAL,       INTENT(IN)  :: condition
#ifndef __PYTHON__
    REAL(sp), DIMENSION(:,:), allocatable :: inverse_sp
#else
    REAL(sp), DIMENSION(size(matrix,2),size(matrix,1)) :: inverse_sp
#endif

#ifndef __PYTHON__
    if (.not. allocated(inverse_sp)) allocate(inverse_sp(size(matrix,2),size(matrix,1)))
#endif
    inverse_sp = real(inverse_dp(real(matrix,dp), condition), sp)

  END FUNCTION inverse_sp

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

  FUNCTION min_diag_dp(matrix, n)

    IMPLICIT NONE

    REAL(dp), DIMENSION(:,:), INTENT(IN) :: matrix
    INTEGER(i4),              INTENT(IN) :: n
#ifndef __PYTHON__
    REAL(dp), DIMENSION(:), allocatable  :: min_diag_dp
#else
    REAL(dp), DIMENSION(size(matrix,1)-abs(n)) :: min_diag_dp
#endif

    if (size(matrix,1) /= size(matrix,2)) stop 'min_diag_dp: array must be squared matrix.'
    if (abs(n)         >= size(matrix,1)) stop 'min_diag_dp: n is to big. You need to choose a minor diagonal within the matrix.'
#ifndef __PYTHON__
    if (.not. allocated(min_diag_dp)) allocate(min_diag_dp(size(matrix,1) - abs(n)))
#endif

    if (n >= 0_i4) then
        min_diag_dp = diag_dp(matrix(1_i4:size(matrix,1)-n,1_i4+n:size(matrix,1)))
    else
        min_diag_dp = diag_dp(matrix(1_i4-n:size(matrix,1),1_i4:size(matrix,1)+n))
    end if

  END FUNCTION min_diag_dp

  FUNCTION min_diag_sp(matrix, n)

    IMPLICIT NONE

    REAL(sp), DIMENSION(:,:), INTENT(IN) :: matrix
    INTEGER(i4),              INTENT(IN) :: n
#ifndef __PYTHON__
    REAL(sp), DIMENSION(:), allocatable  :: min_diag_sp
#else
    REAL(sp), DIMENSION(size(matrix,1)-abs(n)) :: min_diag_sp
#endif

    if (size(matrix,1) /= size(matrix,2)) stop 'min_diag_sp: array must be squared matrix.'
    if (abs(n)         >= size(matrix,1)) stop 'min_diag_sp: n is to big. You need to choose a minor diagonal within the matrix.'
#ifndef __PYTHON__
    if (.not. allocated(min_diag_sp)) allocate(min_diag_sp(size(matrix,1) - abs(n)))
#endif

    if (n >= 0_i4) then
        min_diag_sp = diag_sp(matrix(1_i4:size(matrix,1)-n,1_i4+n:size(matrix,1)))
    else
        min_diag_sp = diag_sp(matrix(1_i4-n:size(matrix,1),1_i4:size(matrix,1)+n))
    end if

  END FUNCTION min_diag_sp

  FUNCTION min_diag_i4(matrix, n)

    IMPLICIT NONE

    INTEGER(i4), DIMENSION(:,:), INTENT(IN) :: matrix
    INTEGER(i4),                 INTENT(IN) :: n
#ifndef __PYTHON__
    INTEGER(i4), DIMENSION(:), allocatable  :: min_diag_i4
#else
    INTEGER(i4), DIMENSION(size(matrix,1)-abs(n)) :: min_diag_i4
#endif

    if (size(matrix,1) /= size(matrix,2)) stop 'min_diag_i4: array must be squared matrix.'
    if (abs(n)         >= size(matrix,1)) stop 'min_diag_i4: n is to big. You need to choose a minor diagonal within the matrix.'
#ifndef __PYTHON__
    if (.not. allocated(min_diag_i4)) allocate(min_diag_i4(size(matrix,1) - abs(n)))
#endif

    if (n >= 0_i4) then
        min_diag_i4 = diag_i4(matrix(1_i4:size(matrix,1)-n,1_i4+n:size(matrix,1)))
    else
        min_diag_i4 = diag_i4(matrix(1_i4-n:size(matrix,1),1_i4:size(matrix,1)+n))
    end if

  END FUNCTION min_diag_i4

  FUNCTION min_diag_i8(matrix, n)

    IMPLICIT NONE

    INTEGER(i8), DIMENSION(:,:), INTENT(IN) :: matrix
    INTEGER(i4),                 INTENT(IN) :: n
#ifndef __PYTHON__
    INTEGER(i8), DIMENSION(:), allocatable  :: min_diag_i8
#else
    INTEGER(i8), DIMENSION(size(matrix,1)-abs(n)) :: min_diag_i8
#endif

    if (size(matrix,1) /= size(matrix,2)) stop 'min_diag_i8: array must be squared matrix.'
    if (abs(n)         >= size(matrix,1)) stop 'min_diag_i8: n is to big. You need to choose a minor diagonal within the matrix.'
#ifndef __PYTHON__
    if (.not. allocated(min_diag_i8)) allocate(min_diag_i8(size(matrix,1) - abs(n)))
#endif

    if (n >= 0_i4) then
        min_diag_i8 = diag_i8(matrix(1_i4:size(matrix,1)-n,1_i4+n:size(matrix,1)))
    else
        min_diag_i8 = diag_i8(matrix(1_i4-n:size(matrix,1),1_i4:size(matrix,1)+n))
    end if

  END FUNCTION min_diag_i8

  FUNCTION min_diag_lgt(matrix, n)

    IMPLICIT NONE

    LOGICAL, DIMENSION(:,:), INTENT(IN) :: matrix
    INTEGER(i4),             INTENT(IN) :: n
#ifndef __PYTHON__
    LOGICAL, DIMENSION(:), allocatable  :: min_diag_lgt
#else
    LOGICAL, DIMENSION(size(matrix,1)-abs(n)) :: min_diag_lgt
#endif

    if (size(matrix,1) /= size(matrix,2)) stop 'min_diag_lgt: array must be squared matrix.'
    if (abs(n)         >= size(matrix,1)) stop 'min_diag_lgt: n is to big. You need to choose a minor diagonal within the matrix.'
#ifndef __PYTHON__
    if (.not. allocated(min_diag_lgt)) allocate(min_diag_lgt(size(matrix,1) - abs(n)))
#endif

    if (n >= 0_i4) then
        min_diag_lgt = diag_lgt(matrix(1_i4:size(matrix,1)-n,1_i4+n:size(matrix,1)))
    else
        min_diag_lgt = diag_lgt(matrix(1_i4-n:size(matrix,1),1_i4:size(matrix,1)+n))
    end if

  END FUNCTION min_diag_lgt

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

  FUNCTION solve_linear_equations_1_dp(lhs, rhs, condition)

    IMPLICIT NONE

    REAL(dp), DIMENSION(:,:), INTENT(IN) :: lhs
    REAL(dp), DIMENSION(:),   INTENT(IN) :: rhs
    LOGICAL,  OPTIONAL,       INTENT(IN) :: condition
#ifndef __PYTHON__
    REAL(dp), DIMENSION(:), allocatable  :: solve_linear_equations_1_dp
#else
    REAL(dp), DIMENSION(size(lhs,1)) :: solve_linear_equations_1_dp
#endif

    INTEGER(i4) :: ii, nZeilen
    real(dp),    dimension(:,:), allocatable :: ilhs          ! internal lhs
    real(dp),    dimension(:),   allocatable :: irhs          ! internal rhs
    real(dp),    dimension(:),   allocatable :: scale_cols    ! scale matrix for better conditioning
    real(dp),    dimension(:),   allocatable :: scale_rows    !              "
    integer(i4), dimension(:),   allocatable :: ipiv          ! needed for dgesv lapack routine
    integer(i4)                              :: info          !              "
    logical :: icondition

    external :: dgesv  ! Lapack routine to compute solution of real system of linear equations

    nZeilen = size(lhs,1)
    if (size(lhs,2) /= nZeilen) stop 'solve_linear_equations_1_dp: left hand side must be squared matrix.'
    if (size(rhs,1) /= nZeilen) stop 'solve_linear_equations_1_dp: right hand side must have same size as left hand side.'
    ! internal arrays
    allocate(ilhs(nZeilen,nZeilen), irhs(nZeilen))
    ilhs = lhs
    irhs = rhs

    if (present(condition)) then
       icondition = condition
    else
       icondition = .true.
    endif

    if (icondition) then
       ! Condition of matrix
       allocate(scale_cols(nZeilen), scale_rows(nZeilen))
       ! Condition columns
       scale_cols(:) = maxval(abs(ilhs(:,:)),1)
       where (scale_cols(:) < tiny(1.0_dp)) scale_cols(:) = 1.0_dp
       scale_cols(:) = 1.0_dp / scale_cols(:)
       forall(ii=1:nZeilen) ilhs(:,ii) = ilhs(:,ii) * scale_cols(ii)
       ! Condition rows
       scale_rows(:) = maxval(abs(ilhs(:,:)),2)
       where (scale_rows(:) < tiny(1.0_dp)) scale_rows(:) = 1.0_dp
       scale_rows(:) = 1.0_dp / scale_rows(:)
       forall(ii=1:nZeilen) ilhs(ii,:) = ilhs(ii,:) * scale_rows(ii)
       irhs(:) = irhs(:) * scale_rows(:)
    endif

    ! solve linear system of equations
    allocate(ipiv(nZeilen))
    call dgesv(nZeilen, 1, ilhs, nZeilen, ipiv, irhs, nZeilen, info)
    if (info /= 0) stop 'solve_linear_equations_1_dp: Solving of linear system did not work.'

#ifndef __PYTHON__
    if (.not. allocated(solve_linear_equations_1_dp)) allocate(solve_linear_equations_1_dp(nZeilen))
#endif
    solve_linear_equations_1_dp = irhs

    if (icondition) then
       ! rescale result
       solve_linear_equations_1_dp(:) = solve_linear_equations_1_dp(:) * scale_cols(:)

       deallocate(scale_cols, scale_rows)
    endif

    deallocate(ilhs, irhs)
    deallocate(ipiv)

  END FUNCTION solve_linear_equations_1_dp

  FUNCTION solve_linear_equations_1_sp(lhs, rhs, condition)

    IMPLICIT NONE

    REAL(sp), DIMENSION(:,:), INTENT(IN) :: lhs
    REAL(sp), DIMENSION(:),   INTENT(IN) :: rhs
    LOGICAL,  OPTIONAL,       INTENT(IN) :: condition
#ifndef __PYTHON__
    REAL(sp), DIMENSION(:), allocatable  :: solve_linear_equations_1_sp
#else
    REAL(sp), DIMENSION(size(rhs,1)) :: solve_linear_equations_1_sp
#endif

#ifndef __PYTHON__
    if (.not. allocated(solve_linear_equations_1_sp)) allocate(solve_linear_equations_1_sp(size(rhs,1)))
#endif
    solve_linear_equations_1_sp = real(solve_linear_equations_1_dp(real(lhs,dp), real(rhs,dp), condition), sp)

  END FUNCTION solve_linear_equations_1_sp

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

  FUNCTION solve_linear_equations_svd_1_dp(lhs, rhs, condition)

    IMPLICIT NONE

    REAL(dp), DIMENSION(:,:), INTENT(IN) :: lhs
    REAL(dp), DIMENSION(:),   INTENT(IN) :: rhs
    LOGICAL,  OPTIONAL,       INTENT(IN) :: condition
#ifndef __PYTHON__
    REAL(dp), DIMENSION(:), allocatable  :: solve_linear_equations_svd_1_dp
#else
    REAL(dp), DIMENSION(size(lhs,1)) :: solve_linear_equations_svd_1_dp
#endif

    INTEGER(i4) :: ii, nZeilen
    real(dp),    dimension(:,:), allocatable :: ilhs          ! internal lhs
    real(dp),    dimension(:),   allocatable :: irhs          ! internal rhs
    real(dp),    dimension(:),   allocatable :: scale_cols    ! scale matrix for better conditioning
    real(dp),    dimension(:),   allocatable :: scale_rows    !              "
    real(dp),    dimension(:),   allocatable :: work          ! needed for dgesvd lapack routine
    real(dp),    dimension(:),   allocatable :: svdw          !              "
    real(dp),    dimension(:,:), allocatable :: svdu          !              "
    real(dp),    dimension(:,:), allocatable :: svdv          !              "
    integer(i4)                              :: lwork         !              "
    integer(i4)                              :: info          !              "
    real(dp),    parameter                   :: svdtol = 1.0e-5_dp ! if <svdtol*maxval(svdw), then set svdw=0
    logical :: icondition

    external :: dgesvd ! Lapack routine to compute singular value decomposition of matrix

    nZeilen = size(lhs,1)
    if (size(lhs,2) /= nZeilen) stop 'solve_linear_equations_svd_1_dp: left hand side must be squared matrix.'
    if (size(rhs,1) /= nZeilen) stop 'solve_linear_equations_svd_1_dp: right hand side must have same size as left hand side.'
    ! internal arrays
    allocate(ilhs(nZeilen,nZeilen), irhs(nZeilen))
    ilhs = lhs
    irhs = rhs

    if (present(condition)) then
       icondition = condition
    else
       icondition = .true.
    endif

    if (icondition) then
       ! Condition of matrix
       allocate(scale_cols(nZeilen), scale_rows(nZeilen))
       ! Condition columns
       scale_cols(:) = maxval(abs(ilhs(:,:)),1)
       where (scale_cols(:) < tiny(1.0_dp)) scale_cols(:) = 1.0_dp
       scale_cols(:) = 1.0_dp / scale_cols(:)
       forall(ii=1:nZeilen) ilhs(:,ii) = ilhs(:,ii) * scale_cols(ii)
       ! Condition rows
       scale_rows(:) = maxval(abs(ilhs(:,:)),2)
       where (scale_rows(:) < tiny(1.0_dp)) scale_rows(:) = 1.0_dp
       scale_rows(:) = 1.0_dp / scale_rows(:)
       forall(ii=1:nZeilen) ilhs(ii,:) = ilhs(ii,:) * scale_rows(ii)
       irhs(:) = irhs(:) * scale_rows(:)
    endif

    ! solve linear system of equations
    ! Use Lapack singular value decomposition and Numerical Recipes solution
    allocate(work(1))
    allocate(svdu(nZeilen,nZeilen), svdv(nZeilen,nZeilen), svdw(nZeilen))
    call dgesvd('A', 'A', nZeilen, nZeilen, ilhs, nZeilen, svdw, svdu, nZeilen, svdv, nZeilen, work, -1, info)
    lwork = int(work(1),i4)
    deallocate(work)
    allocate(work(lwork))
    call dgesvd('A', 'A', nZeilen, nZeilen, ilhs, nZeilen, svdw, svdu, nZeilen, svdv, nZeilen, work, lwork, info)
    if (info /= 0) stop 'solve_linear_equations_svd_1_dp: Solving of linear system did not work.'
    ! write(*,*) 'Matrix condition ', maxval(abs(svdw))/minval(abs(svdw))
    where (svdw < svdtol*maxval(svdw)) svdw = 0.0_dp

#ifndef __PYTHON__
    if (.not. allocated(solve_linear_equations_svd_1_dp)) allocate(solve_linear_equations_svd_1_dp(nZeilen))
#endif
    ! svdv is V**T from dgesvd, svdksb wants V
    solve_linear_equations_svd_1_dp = svdksb(svdu,svdw,transpose(svdv),irhs)

    if (icondition) then
       ! rescale result
       solve_linear_equations_svd_1_dp(:) = solve_linear_equations_svd_1_dp(:) * scale_cols(:)

       deallocate(scale_cols, scale_rows)
    endif

    deallocate(ilhs, irhs)
    deallocate(work, svdu, svdv, svdw)

  END FUNCTION solve_linear_equations_svd_1_dp

  FUNCTION solve_linear_equations_svd_1_sp(lhs, rhs, condition)

    IMPLICIT NONE

    REAL(sp), DIMENSION(:,:), INTENT(IN) :: lhs
    REAL(sp), DIMENSION(:),   INTENT(IN) :: rhs
    LOGICAL,  OPTIONAL,       INTENT(IN) :: condition
#ifndef __PYTHON__
    REAL(sp), DIMENSION(:), allocatable  :: solve_linear_equations_svd_1_sp
#else
    REAL(sp), DIMENSION(size(rhs,1)) :: solve_linear_equations_svd_1_sp
#endif

#ifndef __PYTHON__
    if (.not. allocated(solve_linear_equations_svd_1_sp)) allocate(solve_linear_equations_svd_1_sp(size(rhs,1)))
#endif
    solve_linear_equations_svd_1_sp = real(solve_linear_equations_svd_1_dp(real(lhs,dp), real(rhs,dp), condition), sp)

  END FUNCTION solve_linear_equations_svd_1_sp

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

  FUNCTION solve_linear_equations_band_1_dp(lhsb, rhs, l, u, condition)

    IMPLICIT NONE

    REAL(dp), DIMENSION(:,:), INTENT(IN) :: lhsb !banded form of the lhs
    REAL(dp), DIMENSION(:),   INTENT(IN) :: rhs
    INTEGER(i4),              INTENT(IN) :: l, u
    LOGICAL,  OPTIONAL,       INTENT(IN) :: condition
#ifndef __PYTHON__
    REAL(dp), DIMENSION(:), allocatable  :: solve_linear_equations_band_1_dp
#else
    REAL(dp), DIMENSION(size(rhs,1)) :: solve_linear_equations_band_1_dp
#endif

    INTEGER(i4) :: ii, nSpalten, nZeilen
    real(dp),    dimension(:,:), allocatable :: ilhsb         ! internal lhsb
    real(dp),    dimension(:),   allocatable :: irhs          ! internal rhs
    real(dp),    dimension(:),   allocatable :: scale_cols    ! scale matrix for better conditioning
!    real(dp),    dimension(:),   allocatable :: scale_rows    !              "
    integer(i4), dimension(:),   allocatable :: ipiv          ! needed for dgesv lapack routine
    integer(i4)                              :: info          !              "
    logical :: icondition

    external :: dgbsv  ! Lapack routine to compute solution of real system of linear equations with banded matrix

    nZeilen  = size(lhsb,1)
    nSpalten = size(lhsb,2)
    if (l+u+1_i4    /= nZeilen ) stop 'solve_linear_equations_1_dp: Given banded matrix must have l+u+1 rows.'
    if (size(rhs,1) /= nSpalten) stop 'solve_linear_equations_1_dp: right hand side must have same size as left hand side.'
    ! internal arrays
    allocate(ilhsb(nZeilen + l,nSpalten), irhs(nSpalten))
    ilhsb = 0.0_dp

    ilhsb(1+l:nZeilen + l,:) = lhsb
    irhs                     = rhs

    if (present(condition)) then
       icondition = condition
    else
       icondition = .true.
    endif

    if (icondition) then
       ! Condition of matrix
       allocate(scale_cols(nSpalten))
       ! Condition columns
       scale_cols(:) = maxval(abs(ilhsb(:,:)),1)
       where (scale_cols(:) < tiny(1.0_dp)) scale_cols(:) = 1.0_dp
       scale_cols(:) = 1.0_dp / scale_cols(:)
       forall(ii=1:nSpalten) ilhsb(:,ii) = ilhsb(:,ii) * scale_cols(ii)
    endif

    ! solve linear system of equations
    allocate(ipiv(nSpalten))
    call dgbsv(nSpalten, l, u, 1, ilhsb, nZeilen+l, ipiv, irhs, nSpalten, info)
    if (info /= 0) stop 'solve_linear_equations_band_1_dp: Solving of linear system did not work.'

#ifndef __PYTHON__
    if (.not. allocated(solve_linear_equations_band_1_dp)) allocate(solve_linear_equations_band_1_dp(nSpalten))
#endif
    solve_linear_equations_band_1_dp = irhs

    if (icondition) then
       ! rescale result
       solve_linear_equations_band_1_dp(:) = solve_linear_equations_band_1_dp(:) * scale_cols(:)

       deallocate(scale_cols)
    endif

    deallocate(ilhsb, irhs)
    deallocate(ipiv)

  END FUNCTION solve_linear_equations_band_1_dp

  FUNCTION solve_linear_equations_band_1_sp(lhsb, rhs, l, u, condition)

    IMPLICIT NONE

    REAL(sp), DIMENSION(:,:), INTENT(IN) :: lhsb !banded form of the lhs
    REAL(sp), DIMENSION(:),   INTENT(IN) :: rhs
    INTEGER(i4),              INTENT(IN) :: l, u
    LOGICAL,  OPTIONAL,       INTENT(IN) :: condition
#ifndef __PYTHON__
    REAL(sp), DIMENSION(:), allocatable  :: solve_linear_equations_band_1_sp
#else
    REAL(sp), DIMENSION(size(rhs,1)) :: solve_linear_equations_band_1_sp
#endif

#ifndef __PYTHON__
    if (.not. allocated(solve_linear_equations_band_1_sp)) allocate(solve_linear_equations_band_1_sp(size(rhs,1)))
#endif
    solve_linear_equations_band_1_sp = real(solve_linear_equations_band_1_dp(real(lhsb,dp), real(rhs,dp), l, u, condition), sp)

  END FUNCTION solve_linear_equations_band_1_sp

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

  ! From Numerical Recipes in F77 book
  ! Solves A*X=B for a vector X, where A is specified by the arrays u, w, v as returned by
  ! svdcmp. m and n are the logical dimensions of a, and will be equal for square matrices. mp
  ! and np are the physical dimensions of a. b(1:m) is the input right-handside. x(1:n) is
  ! the output solution vector. No input quantities are destroyed, so the routine may be called
  ! sequentially with different bs.
  FUNCTION svdksb_dp(u,w,v,b)

    IMPLICIT NONE

    REAL(dp), DIMENSION(:,:), INTENT(IN) :: u,v
    REAL(dp), DIMENSION(:),   INTENT(IN) :: w,b
    REAL(dp), DIMENSION(size(u,2))       :: svdksb_dp

    REAL(dp), DIMENSION(size(u,2)) :: tmp
    REAL(dp), DIMENSION(size(w))   :: wtmp

    wtmp = w
    where (abs(w) < tiny(1.0_dp)) wtmp = 1.0_dp
    tmp = matmul(b,u) / wtmp
    where (abs(w) < tiny(1.0_dp)) tmp = 0.0_dp
    svdksb_dp = matmul(v,tmp)

  END FUNCTION svdksb_dp

  FUNCTION svdksb_sp(u,w,v,b)

    IMPLICIT NONE

    REAL(sp), DIMENSION(:,:), INTENT(IN) :: u,v
    REAL(sp), DIMENSION(:),   INTENT(IN) :: w,b
    REAL(sp), DIMENSION(size(u,2))       :: svdksb_sp

    REAL(sp), DIMENSION(size(u,2)) :: tmp
    REAL(sp), DIMENSION(size(w))   :: wtmp

    wtmp = w
    where (abs(w) < tiny(1.0_sp)) wtmp = 1.0_sp
    tmp = matmul(b,u) / wtmp
    where (abs(w) < tiny(1.0_sp)) tmp = 0.0_sp
    svdksb_sp = matmul(v,tmp)

  END FUNCTION svdksb_sp

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

END MODULE mo_linear_algebra


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE SPTXYZ (XN, YN, XN3, YN3, ZN3, IXNP, NRNP, ZCORD,
     &     NSPL, ZS, XS, XS2, YS, YS2, SCR )
C=======================================================================

C     --*** SPLTYZ *** (GEN3D) Calculate 3D coordinates for Spline translation
C     --   Written by Greg Sjaardema - 02/06/89
C     --
C     --Parameters:
C     --   XN, YN - IN - the 2D coordinates, destroyed
C     --   XN3, YN3, ZN3 - OUT - the 3D coordinates
C     --   IXNP - IN - the new index for each node
C     --   NRNP - IN - the number of new nodes generated for each node
C     --   ZCORD - SCRATCH - size = NNREPL, holds z coordinate for transformations
C     --   NSPL - IN - number of spline points
C     --   ZS, XS, YS - IN - spline data points
C     --   XS2, YS2 - SCRATCH - temporary storage for spline data
C     --   SCR - SCRATCH - temporary storage for spoine data
C     --
C     --Common Variables:
C     --   Uses NDIM, NUMNP of /DBNUMS/
C     --   Uses NDIM3, NUMNP3 of /DBNUM3/
C     --   Uses DOTRAN, NNREPL, DIM3, NRTRAN, D3TRAN, ZGRAD,
C     --      CENTER, NUMCOL, NUMROW of /PARAMS/

      INCLUDE 'g3_dbnums.blk'
      INCLUDE 'g3_dbnum3.blk'
      INCLUDE 'g3_params.blk'

      PARAMETER (BINGO = 1.0E38)
      PARAMETER (TOLER = 1.0E-8)

      REAL XN(NUMNP), YN(NUMNP),
     &     XN3(NUMNP3), YN3(NUMNP3), ZN3(NUMNP3)
      INTEGER IXNP(*), NRNP(*)
      REAL ZCORD(NNREPL)
      REAL ZS(NSPL), XS(NSPL), XS2(NSPL), YS(NSPL), YS2(NSPL)
      REAL SCR(NSPL), SLLFT(2), SLRGT(2)

C     ... CALCULATE THE THICKNESS INCREMENT FOR EACH TRANSLATION

      IBLK = 0
      ZTOT = 0.0
   10 CONTINUE
      IBLK = IBLK + 1
      IF (NRTRAN(IBLK) .GT. 0) THEN
         ZTOT = ZTOT + D3TRAN(IBLK)
         IF (IBLK .LT. MAXINT) GO TO 10
      END IF

      NXTNR = 1
      IBLK = 0
      ZEND = 0.0

   20 CONTINUE
      IBLK = IBLK + 1
      IF (NRTRAN(IBLK) .GT. 0) THEN
         ZBEG = ZEND
         ZEND = ZBEG + D3TRAN(IBLK)
         CALL INIGRD (ZBEG/ZTOT, ZEND/ZTOT, ZGRAD(IBLK),
     *        NRTRAN(IBLK), NRTRAN(IBLK)+1, ZCORD(NXTNR) )
         NXTNR = NXTNR + NRTRAN(IBLK)
         IF (IBLK .LT. MAXINT) GO TO 20
      END IF

      CALL SPLINE (ZS, XS, NSPL, SLLFT(1), SLRGT(1), XS2, SCR)
      CALL SPLINE (ZS, YS, NSPL, SLLFT(2), SLRGT(2), YS2, SCR)

      ZMAX = ZS(NSPL)

      KLO = 1

      DO 80 NR = 1, NNREPL
         Z = ZCORD(NR) * ZMAX
         CALL HUNT (ZS, NSPL, Z, KLO)
         KLO = MIN(NSPL-1, MAX(1, KLO))

         H  = ZS(KLO+1) - ZS(KLO)
         A  = (ZS(KLO+1)-Z)/H
         B  = (Z-ZS(KLO))/H

         XOFF = A * XS(KLO) + B * XS(KLO+1) +
     *        ((A**3-A) * XS2(KLO)+(B**3-B)*XS2(KLO+1)) * (H**2) / 6.

         YOFF = A * YS(KLO) + B * YS(KLO+1) +
     *        ((A**3-A) * YS2(KLO)+(B**3-B)*YS2(KLO+1)) * (H**2) / 6.

         DO 70 INP = 1, NUMNP
            JNP = IXNP(INP) - 1
            XN3(JNP+NR) = XN(INP) + XOFF
            YN3(JNP+NR) = YN(INP) + YOFF
            ZN3(JNP+NR) = Z
   70    CONTINUE
   80 CONTINUE

      RETURN
      END

max_stars_repo_name: tokusanya/seacas

max_issues_repo_name: tokusanya/seacas

max_forks_repo_name: tokusanya/seacas


text:
c *****************************************************
c    diskmetg.f
c 26 Jan., 2011   written by D. Kawata
c *****************************************************
c70*******************************************************************

      subroutine diskmetg(ngr,ngz,dfeh,feh0,dalfe,alfe0
     &  ,nhhalo,fehhalo)
      include './define.f'
      integer ngr,ngz
      double precision dfeh,feh0,dalfe,alfe0
      integer i,j,pn
      double precision rp,fehp,alfep,lmetp,metp
      double precision xsolh,xsolz,xh0,fehhalo,nhp,nhhalo

      xh0=0.76d0
      xsolh=0.706d0
      xsolz=0.019d0

      if(myrank.eq.0) then
        write(6,*) 'nhhalo,fehhalo=',nhhalo,fehhalo
      endif

      do i=0,ngr
        do j=0,ngz
          pn=id_rzg(i,j)
c check density
          fh_rzg(pn)=(xh0+(xsolh-xh0)*(10.0d0**fehhalo))
          nhp=rho_rzg(pn)*fh_rzg(pn)*(DU/MP)
          if(nhp.lt.nhhalo) then
            fehp=fehhalo
          else
c *** radius
            rp=x_rzg(pn)
            fehp=dfeh*(rp*LUKPC)+feh0
          endif
          alfep=dalfe*fehp+alfe0
c *** set [Z/H] from [Fe/H] and [al/Fe]
          lmetp=fehp
          if(fehp.lt.fehp+alfep) then
c if alpha overabundant
            lmetp=fehp+alfep
          endif
          if(fehp.lt.fehp-alfep) then
c for the case of C or N enhanced
            lmetp=fehp-alfep
          endif
          metp=10.0d0**lmetp
          fh_rzg(pn)=(xh0+(xsolh-xh0)*(10.0d0**fehp))
          lmet_rzg(pn)=dlog10((1.0d0/xsolh)
     &        *(xh0+(xsolh-xh0)*metp)*metp)
        enddo
      enddo

      return
      end

max_stars_repo_name: kawatadaisuke/diskm-buildini

max_issues_repo_name: kawatadaisuke/diskm-buildini

max_forks_repo_name: kawatadaisuke/diskm-buildini


text:
c-----DMFIT, v 0.1
c
c     this test program tests the dmfit subroutine and function
c     it also shows an application to KK dark matter

      program test

      implicit none

      character*32 arg
c-----variables needed to run the routines: mass
      real*8 MX,z
c     annihilation channel, smoothing parameter, switch diff/integral (see below)
      integer CH,HSM,SWITCH,ii
c     remember that the following channels, corresponding to WIMP annihilation modes
c     are implemented:
c
c     1: e^+e^-
c     2: mu^+ mu^-
c     3: tau^+ tau^-
c     4: b \bar b
c     5: t \bar t
c     6: gluon gluon
c     7: W+ W-
c     8: Z Z
c     9: c \bar c

c     gamma ray energy, in GeV
      real*8 EE

c-----outputs:
      real*8 DNDE,DNDM

c-----the actual dmfit functions
      real*8 dmfit_de,dmfit_dm

c-----test case: a 100 GeV WIMP in b\bar b, same smoothing as DarkSUSY's default
      CALL getarg(1,arg)
      read(arg,*) MX
      CALL getarg(2,arg)
      read(arg,*) CH
c      MX=6.d0
c      CH=12
      HSM=0

c-----let's look at an energy of 10 GeV
      EE=0.2d0

      call dmfit_load ("new_gammamc_dif.dat")

c     now let's get the integrated flux

c$$$      DNDE=dmfit_de(MX,CH,EE)
c$$$      DNDM=dmfit_dm(MX,CH,EE)

      write(*,*) '# mass',MX
      write(*,*) '# channel',CH
c      write(*,*) 'energy',EE

c$$$      write(*,*)
c$$$      write(*,*) 'dN/dE:',DNDE
c$$$      write(*,*) 'd(dN/dE)/dM:',DNDM
      write(*,*) '# idx energy Log(E/MX) EdN/dE d(dN/dE)/dM'
      do ii=0,249
         z=(ii+0.5)/250.
         EE=MX*10.**(10*(z-1.))
         DNDE=dmfit_de(MX,CH,EE)
c          write(*,*) DNDE
c         DNDM=dmfit_dm(MX,CH,EE)
         write(*,*) ii, EE,log10(EE)-log10(MX),DNDE!, DNDM
      enddo

      stop
      end

c-----these are the routines needed to be included:

      include 'dmfit_load.f'
      include 'dmfit_func.f'
      include 'dmfit_comm.f'

max_stars_repo_name: fermi-lat/Likelihood

max_issues_repo_name: fermi-lat/Likelihood

max_forks_repo_name: fermi-lat/Likelihood


text:
	SUBROUTINE IS_CNTM ( timstr, iotarr, strtim, iret )
C************************************************************************
C* IS_CNTM 								*
C*									*
C* This subroutine decodes the cancellation time string.		*
C*                                                                      *
C* IS_CNTM ( TIMSTR, IOTARR, STRTIM, IRET )				*
C*									*
C* Input parameters:							*
C*	TIMSTR		CHAR*		Bulletin time string		*
C*	IOTARR (5)	INTEGER		Bull. time - YYYY,MM,DD,HH,MM   *
C*									*
C* Output parameters:							*
C*	STRTIM		CHAR*		GEMPAK time			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -2 = format error		*
C*									*
C**									*
C* Log:									*
C* F. J. Yen/NCEP	 2/02	Created from part of IS_PRMB.		*
C************************************************************************
C*
	CHARACTER*(*)	timstr, strtim
	INTEGER  	iotarr (*)
C*
	INTEGER  	itarr (5)
C------------------------------------------------------------------------
	iret  = 0
	CALL ST_LSTR ( timstr, lens, ier )
	IF ( ( lens .eq. 5 ) .or. ( lens .eq. 7 ) ) THEN

	    IF ( timstr ( lens:lens ) .eq. 'Z' .or.
     +		    timstr ( lens:lens ) .eq. '=' .or.
     +		    timstr ( lens:lens ) .eq. '.' )
     +			lens = lens - 1
	  ELSE IF ( lens .eq. 8 ) THEN
	    IF ( timstr ( lens - 1:lens ) .eq. 'Z.' )
     +			lens = lens - 2
	END IF
C
	IF ( ( lens .eq. 4 ) .or. ( lens .eq. 6 ) ) THEN
	    DO i = 1, 3
	        itarr ( i ) = iotarr ( i )
	    END DO
	    CALL ST_INTG ( timstr ( lens - 3:lens ),
     +			   ihhmm, ier )
	    IF ( ier .lt. 0 ) THEN
		iret = -2
		RETURN
	    END IF
C
	    itarr ( 4 ) = ihhmm / 100
	    itarr ( 5 ) = MOD ( ihhmm, 100 )
	    CALL TI_MDIF ( itarr, iotarr, nmin, ier )
	    IF ( nmin .lt. ( -10 ) ) THEN
		CALL TI_ADDD ( itarr, itarr, ier )
	      ELSE IF ( nmin .gt. 1380 ) THEN
		CALL TI_SUBD ( itarr, itarr, ier )
	    END IF
	    IF ( ier .lt. 0 ) THEN
		iret = -2
		RETURN
	    END IF
C
	    CALL TI_ITOC ( itarr, strtim, ier )
	    IF ( ier .lt. 0 ) THEN
		iret = -2
		RETURN
	    END IF
	  ELSE
	    iret = -2
	    RETURN
	END IF
C*
	RETURN
	END


text:
Claudia Krich taught in the Davis elementary schools from 1990 until she retired in 2010.  She led the Willett Elementary School Willett Elementary Choir from 2001 during that time, and continues to do so after retirement.  She received the Harmony in Our Lives Award in 2010 for her work in school music.

Source:

Goss, Kelly; 2010 Oct. 20, http://www.davisenterprise.com/ArchivedStories0/Harmony_is_her_life/ Harmony is her life, Davis Enterprise.  Retrieved 25 Jan. 2013.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*###[ ffcxyz:
	subroutine ffcxyz(cy,cz,cdyz,cd2yzz,ivert,sdelpp,sdelps,etalam,
     +			etami,delps,xpi,piDpj,isoort,ldel2s,ns,ier)
***#[*comment:***********************************************************
*									*
*	calculate in a numerically stable way				*
*									*
*	cz(1,2) = (-p(ip1).p(is2) +/- sdelpp)/xpi(ip1)			*
*	cy(1,2) = (-p(ip1).p(is2) +/- sdisc)/xpi(ip1)			*
*			disc = slam1 + 4*eta*xpi(ip)/slam		*
*									*
*	cy(3,4) = 1-cy(1,2)						*
*	cz(3.4) = 1-cz(1,2)						*
*	cdyz(i,j) = cy(i) - cz(j)					*
*									*
*	Input:	ivert		(integer)	1,2 of 3		*
*		sdelpp		(real)		sqrt(lam(p1,p2,p3))/2	*
*		sdelps		(real)		sqrt(-lam(p,ma,mb))/2	*
*		etalam		(real)		det(si.sj)/det(pi.pj)	*
*		etami(6)	(real)		si.si - etalam		*
*		xpi(ns)		(real)		standard		*
*		piDpj(ns,ns)	(real)		standard		*
*		ns		(integer)	dim of xpi,piDpj	*
*									*
*	Output:	cy(4),cz(4),cdyz(4,4)	(complex)	see above	*
*									*
*	Calls:	??							*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments:
*
	integer ivert,isoort(2),ns,ier
	logical ldel2s
	DOUBLE COMPLEX cy(4),cz(4),cdyz(2,2),cd2yzz
	DOUBLE PRECISION sdelpp,sdelps,etalam,etami(6),delps,xpi(ns),
     +		piDpj(ns,ns)
*
*	local variables:
*
	integer i,j,ip1,ip2,ip3,is1,is2,is3
	DOUBLE COMPLEX c
	DOUBLE PRECISION absc,y(4)
	DOUBLE PRECISION delps1,disc,hulp,xlosn
*
*	common blocks:
*
	include 'ff.h'
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations:
*  #[ set up pointers:
	if ( ldel2s .and. ivert .ne. 1 ) goto 100
	is1 = ivert
	is2 = ivert+1
	if ( is2 .eq. 4 ) is2 = 1
	is3 = ivert-1
	if ( is3 .eq. 0 ) is3 = 3
	ip1 = is1 + 3
*	ip2 = is2 + 3
*	ip3 = is3 + 3
	isoort(1) = -10
	isoort(2) = -10
*  #] set up pointers:
*  #[ test input:
	if ( ltest .and. xpi(ip1) .eq. 0 ) then
	    call fferr(47,ier)
	    return
	endif
*  #] test input:
*  #[ get cypm,czpm:
	hulp = sdelps/xpi(ip1)
	cz(1) = DCMPLX(piDpj(ip1,is2)/xpi(ip1),-hulp)
	cz(2) = DCMPLX(piDpj(ip1,is2)/xpi(ip1),+hulp)
	disc = delps/sdelpp
	call ffroot(y(1),y(2),xpi(ip1),piDpj(ip1,is2),etami(is2),disc,
     +								ier)
	cy(1) = y(1)
	cy(2) = y(2)
*  #] get cypm,czpm:
*  #[ get cypm1,czpm1:
	if ( xpi(is1) .eq. xpi(is2) ) then
	    cy(4) = cy(1)
	    cy(3) = cy(2)
	    cz(4) = cz(1)
	    cz(3) = cz(2)
	else
	    cz(3) = 1 - cz(1)
	    cz(4) = 1 - cz(2)
	    if ( absc(cz(3)).lt.xloss .or. absc(cz(4)).lt.xloss ) then
		cz(3) =DCMPLX(-piDpj(ip1,is1)/xpi(ip1),+hulp)
		cz(4) =DCMPLX(-piDpj(ip1,is1)/xpi(ip1),-hulp)
	    endif
	    y(3) = 1 - y(1)
	    y(4) = 1 - y(2)
	    if ( abs(y(3)) .lt. xloss .or. abs(y(4)) .lt. xloss ) then
		call ffroot(y(4),y(3),xpi(ip1),-piDpj(ip1,is1),
     +						etami(is1),disc,ier)
	    endif
	    cy(3) = y(3)
	    cy(4) = y(4)
	endif
*  #] get cypm1,czpm1:
*  #[ get cdypzp, cdypzm:
	cdyz(2,1) = DCMPLX(disc/xpi(ip1),+hulp)
	cdyz(2,2) = DCMPLX(disc/xpi(ip1),-hulp)
	cdyz(1,1) = -cdyz(2,2)
	cdyz(1,2) = -cdyz(2,1)
	cd2yzz = 2*disc/xpi(ip1)
	goto 200
*  #] get cdypzp, cdypzm:
*  #[ special case, get indices:
  100	continue
	if ( ivert.eq.2 ) then
	    is1 = 2
	    ip1 = 5
	else
	    is1 = 1
	    ip1 = 6
	endif
	isoort(1) = -100
	isoort(2) = -100
*  #] special case, get indices:
*  #[ get cypm,czpm:
*
*	special case del2s = 0, hence the roots are not the real roots
*	but z_2'' = (z_2'-1)/delta, z''_3 = -z'_3/delta
*
	hulp = sdelps/xpi(3)
	disc = delps/sdelpp
	if ( ivert .eq. 3 ) then
	    hulp = -hulp
	    disc = -disc
	endif
	cz(1) = DCMPLX(piDpj(is1,3)/xpi(3),-hulp)
	cz(2) = DCMPLX(piDpj(is1,3)/xpi(3),+hulp)
	call ffroot(y(1),y(2),xpi(3),piDpj(is1,3),etami(is1),disc,ier)
	cy(1) = y(1)
	cy(2) = y(2)
*  #] get cypm,czpm:
*  #[ get cypm1,czpm1:
	cz(3) = 1 - cz(1)
	cz(4) = 1 - cz(2)
	if ( absc(cz(3)).lt.xloss .or. absc(cz(4)).lt.xloss ) then
	    if ( lwrite ) print *,'cz(3,4) = ',cz(3),cz(4)
	    if ( ivert.eq.2 ) then
		cz(3) =DCMPLX(piDpj(ip1,3)/xpi(3),+hulp)
		cz(4) =DCMPLX(piDpj(ip1,3)/xpi(3),-hulp)
	    else
		cz(3) =DCMPLX(-piDpj(ip1,3)/xpi(3),+hulp)
		cz(4) =DCMPLX(-piDpj(ip1,3)/xpi(3),-hulp)
	    endif
	    if ( lwrite ) print *,'cz(3,4)+= ',cz(3),cz(4)
	endif
	y(3) = 1 - y(1)
	y(4) = 1 - y(2)
	if ( abs(y(3)) .lt. xloss .or. abs(y(4)) .lt. xloss ) then
	    if ( lwrite ) print *,'y(3,4) = ',y(3),y(4)
	    if ( ivert .eq. 2 ) then
		call ffroot(y(4),y(3),xpi(3),piDpj(ip1,3),etami(ip1),
     +							disc,ier)
	    else
		call ffroot(y(4),y(3),xpi(3),-piDpj(ip1,3),etami(ip1),
     +							disc,ier)
	    endif
	    if ( lwrite ) print *,'y(3,4)+= ',y(3),y(4)
	endif
	cy(3) = y(3)
	cy(4) = y(4)
*  #] get cypm1,czpm1:
*  #[ get cdypzp, cdypzm:
	cdyz(2,1) = DCMPLX(disc/xpi(3),+hulp)
	cdyz(2,2) = DCMPLX(disc/xpi(3),-hulp)
	cdyz(1,1) = -cdyz(2,2)
	cdyz(1,2) = -cdyz(2,1)
	cd2yzz = 2*disc/xpi(3)
*  #] get cdypzp, cdypzm:
*  #[ test output:
  200	continue
	if ( ltest ) then
	    xlosn = xloss**2*DBLE(10)**(-mod(ier,50))
	    do 99 i=1,2
		if ( xlosn*absc(cy(i)+cy(i+2)-1) .gt. precc*max(absc(
     +		    cy(i)),absc(cy(i+2)),x1)) print *,'ffcxyz: error: ',
     +		    'cy(',i+2,')<>1-cy(',i,'):',cy(i+2),cy(i),cy(i+2)+
     +		    cy(i)-1
		if ( xlosn*absc(cz(i)+cz(i+2)-1) .gt. precc*max(absc(
     +		    cz(i)),absc(cz(i+2)),x1)) print *,'ffcxzz: error: ',
     +		    'cz(',i+2,')<>1-cz(',i,'):',cz(i+2),cz(i),cz(i+2)+
     +		    cz(i)-1
		do 98 j=1,2
		    if ( xlosn*absc(cdyz(i,j)-cy(i)+cz(j)) .gt. precc*
     +			max(absc(cdyz(i,j)),absc(cy(i)),absc(cz(j))) )
     +			print *,'ffcxyz: error: cdyz(',i,j,') <> cy(',i,
     +			')-cz(',j,'):',cdyz(i,j),cy(i),cz(j),cdyz(i,j)-
     +			cy(i)+cz(j)
   98		continue
   99	    continue
	endif
*  #] test output:
*###] ffcxyz:
	end
*###[ ffcdwz:
	subroutine ffcdwz(cdwz,cw,cz,i1,j1,l,calpha,calph1,cpi,cdpipj,
     +					cpiDpj,csdeli,csdel2,ns,ier)
***#[*comment:***********************************************************
*									*
*	Recalculate cdwz(i1,j1) = cw(i1) - cz(j1)			*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments:
*
	integer i1,j1,l,ns,ier
	DOUBLE COMPLEX cdwz(2,2),cw(4),cz(4),calpha,calph1,cpi(ns)
	DOUBLE COMPLEX cdpipj(ns,ns),cpiDpj(ns,ns),csdeli(3),csdel2
*
*	local variables:
*
	integer i,n
	DOUBLE COMPLEX cs(8),csum,cfac,c,cddel
	DOUBLE PRECISION xmax,absc,afac
*
*	common blocks:
*
	include 'ff.h'
*
*	statement function
*
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations:
*  #[ calculations:
	if ( l .eq. 1 ) then
	    if ( j1 .eq. 1 ) then
		if ( absc(csdeli(1)+csdel2) .lt. xloss*absc(csdel2) )
     +								then
*		    for example in e-> e g* with eeg loop
*		    first get the difference of csdeli(1) and csdel2:
		    cs(1) = cpi(4)*cdpipj(2,5)
		    cs(2) = -cpiDpj(4,3)*cpiDpj(4,2)
		    cs(3) = cpiDpj(4,3)*cpiDpj(4,5)
		    csum = cs(1)+cs(2)+cs(3)
		    xmax = max(absc(cs(1)),absc(cs(2)),absc(cs(3)))
		    if ( absc(csum) .lt. xloss*xmax ) then
			if ( lwrite ) print *,'ffcdwz: canc in cddel'
			ier = 1
			goto 5
		    endif
		    cddel = csum/(csdel2-csdeli(1))
		    if ( i1 .eq. 1 ) then
			cs(1) = cpi(4)*csdeli(2)
		    else
			cs(1) = -cpi(4)*csdeli(2)
		    endif
		    cs(2) = cddel*cpiDpj(4,2)
		    cs(3) = -cpiDpj(4,3)*csdeli(1)
		    cs(4) = cpiDpj(4,3)*cpiDpj(4,5)
		    cs(5) = -cpi(4)*cpiDpj(5,3)
		    cs(6) = -cddel*csdel2
		    n = 6
		else
		    if ( lwrite ) print *,'ffcdwz: ',
     +			'cannot handle this case yet'
		    ier = ier + 100
		    goto 5
		endif
		csum = 0
		xmax = 0
		do 1 i=1,n
		    csum = csum + cs(i)
		    xmax = max(xmax,absc(cs(i)))
    1		continue
		if ( absc(csum) .lt. xloss*xmax ) then
		    if ( lwrite ) print *,'ffcdwz: still cancellations',
     +			' in cdwz(',i1,j1,l,'): ',csum,xmax
		    ier = ier + 1
		endif
		if (lwrite) print *,'  cdwz(',i1,j1,l,')  =',cdwz(i1,j1)
     +				,min(absc(cw(i1)),absc(cw(i1+2)))
		cdwz(i1,j1) = csum/calph1/cpi(4)/cpi(5)
		if ( cdwz(i1,j1) .eq. 0 .and. csum .ne. 0 ) then
		    print *,'?#$&!! cdwz = 0 but csum != 0, try again'
		    afac = 1/absc(csum)
		    csum = csum*DBLE(afac)
		    cdwz(i1,j1) = csum/calph1/cpi(4)/cpi(5)
		    afac = 1/afac
		    cdwz(i1,j1) = cdwz(i1,j1)*DBLE(afac)
		endif
		if (lwrite) print *,'  cdwz(',i1,j1,l,')+ =',cdwz(i1,j1)
     +				,xmax/absc(calph1*cpi(4)*cpi(5))
	    else
		if ( lwrite ) print *,'ffcdwz:  warning: cannot handle',
     +			' this case cdwz(',i1,j1,l,') yet'
		ier = ier + 100
	    endif
    5	    continue
	elseif ( l .eq. 3 ) then
	    if ( (i1.eq.2 .and. j1.eq.1) .or. (i1.eq.1 .and. j1.eq.2 ) )
     +			then
		cfac = 1/(csdeli(2) + csdeli(3))
		cs(1) = cdpipj(6,5)*cz(j1)
		cs(2) = -calph1*cpi(5)*cz(j1+2)
		if ( max(absc(cdpipj(2,1)),absc(cdpipj(5,6))) .lt.
     +		     max(absc(cdpipj(2,6)),absc(cdpipj(5,1))) ) then
		    cs(3) = cdpipj(2,1)/2
		    cs(4) = cdpipj(5,6)/2
		else
		    cs(3) = cdpipj(2,6)/2
		    cs(4) = cdpipj(5,1)/2
		endif
		cs(5) = cpiDpj(4,3)*cpiDpj(5,3)*cfac
		cs(6) = -cpiDpj(4,3)*cpiDpj(6,3)*cfac
		cs(7) = cpi(3)*cdpipj(5,6)*cfac
		if ( i1 .eq. 1 ) then
		    csum = cs(1)+cs(2)+cs(3)+cs(4) - (cs(5)+cs(6)+cs(7))
		else
		    csum = cs(1)+cs(2)+cs(3)+cs(4) + cs(5)+cs(6)+cs(7)
		endif
		xmax = absc(cs(1))
		do 10 i=2,7
		    xmax = max(xmax,absc(cs(i)))
   10		continue
		if ( absc(csum) .lt. xloss*xmax ) then
*		    this result is not used if it is not accurate (see
*		    ffxc0p)
		    if ( lwrite ) then
			call ffwarn(78,ier,absc(csum),xmax)
		    else
			ier = ier + 1
		    endif
		    xmax = xmax/absc(calpha*cpi(5))
		    if ( xmax .lt. min(absc(cz(j1)),absc(cz(j1+2))) )
     +								then
			if (lwrite) print *,'  cdwz(',i1,j1,l,')  = ',
     +			    cdwz(i1,j1),min(absc(cw(i1)),absc(cw(i1+2)))
			cdwz(i1,j1) = csum/(calpha*cpi(5))
			if (lwrite) print *,'  cdwz(',i1,j1,l,')+ = ',
     +						cdwz(i1,j1),xmax
		    endif
		else
		    if (lwrite) print *,'  cdwz(',i1,j1,l,')  = ',
     +			    cdwz(i1,j1),min(absc(cw(i1)),absc(cw(i1+2)))
		    cdwz(i1,j1) = csum/(calpha*cpi(5))
		    if (lwrite) print *,'  cdwz(',i1,j1,l,')+ = ',
     +				cdwz(i1,j1),xmax/absc(calpha*cpi(5))
		endif
	    else
		if ( lwrite ) print *,'ffcdwz:  warning: cannot handle',
     +			' this case cdwz(',i1,j1,l,') yet'
		ier = ier + 100
	    endif
	else
	    if ( lwrite ) print *,'ffcdwz:  error: l <> 1 or 3 but ',l
	    ier = ier + 100
	endif
*  #] calculations:
*  #[ test output:
	if ( ltest .and. ier .eq. 0 ) then
	    if ( xloss**2*absc(cdwz(i1,j1)-cw(i1)+cz(j1)) .gt. precc*
     +		max(absc(cdwz(i1,j1)),absc(cw(i1)),absc(cz(j1))) )
     +		print *,'ffcdwz:  error: cdwz(',i1,j1,l,') <> cw - cz :'
     +		,cdwz(i1,j1),cw(i1),cz(j1),cw(i1)-cz(j1),
     +		cdwz(i1,j1)-cw(i1)+cz(j1)
	    if ( xloss**2*absc(cdwz(i1,j1)+cw(i1+2)-cz(j1+2)) .gt.
     +		precc*max(absc(cdwz(i1,j1)),absc(cw(i1+2)),
     +		absc(cz(j1+2))) ) print *,'ffcdwz:  error: cdwz(',i1,j1,
     +		l,') <> cz1- cw1:',cdwz(i1,j1),cz(i1+2),cw(j1+2),
     +		cz(i1+2)-cw(j1+2),cdwz(i1,j1)+cw(i1+2)-cz(j1+2)
	    endif
*  #] test output:
*###] ffcdwz:
	end


text:
*----------------------------------------------------------------------*
      subroutine set_r12intm(op,name,dagger,
     &     min_rank,max_rank,ncadiff,iformal,op_info,orb_info)
*----------------------------------------------------------------------*
*     set up intermediate for R12 with open lines within
*     settings analoguous to R12coeff
*     hpvx_mnmx,irestr are chosen appropriately
*----------------------------------------------------------------------*
      implicit none
      include 'stdunit.h'
      include 'opdim.h'
c      include 'def_operator.h'
c      include 'def_operator_array.h'
      include 'mdef_operator_info.h'
      include 'def_orbinf.h'
      include 'par_opnames_gen.h'
      include 'ifc_input.h'

      integer, parameter ::
     &     ntest = 00

      type(operator), intent(inout) ::
     &     op
      character, intent(in) ::
     &     name*(*)
      logical, intent(in) ::
     &     dagger
      integer, intent(in) ::
     &     min_rank, max_rank, ncadiff, iformal

      type(operator_info), intent(inout) ::
     &     op_info
      type(orbinf), intent(inout) ::
     &     orb_info

      integer ::
     &     idx, idx_dum, i, j
      integer ::
     &     occ0(ngastp,2)
      logical ::
c     &     r12fix
     &     cclog
      type(operator), pointer ::
     &     dum1_pnt, dum2_pnt, dum3_pnt
      type(operator_array), pointer ::
     &     opscr(:)

      character(5) ::
     &     op_dummy1*5 = 'DUM-1',
     &     op_dummy2*5 = 'DUM-2',
     &     op_dummy3*5 = 'DUM-3'

      integer, external ::
     &     idx_oplist2

c      ! Do we want amplitudeless intermediates?
c      call get_argument_value('method','CC',lval=cclog)
      cclog = is_keyword_set('method.CC').gt.0

      if(ntest.ge.100)then
        write(lulog,*) '-------------'
        write(lulog,*) ' set_r12intm '
        write(lulog,*) '-------------'
        write(lulog,*) 'Intermediate: ',trim(name)
c        if(r12fix) write(lulog,*)'Fully contracted intermediate'
      endif

c      if(.not.r12fix)then
        allocate(opscr(2))
        allocate(opscr(1)%op)
        allocate(opscr(2)%op)

        ! no external lines, so:
        occ0 = 0
        call set_uop(opscr(1)%op,'scr1',.false.,
     &       occ0,1,orb_info)

        ! the internal lines are defined by a coefficient-like
        ! operator, so we use ...
        call set_r12c(opscr(2)%op,'scr2',dagger,
     &       min_rank,max_rank,ncadiff,iformal,orb_info)

        ! ... and define the intermediate:
        call set_gen_intermediate(op,name,
     &       opscr,2,orb_info)

        call dealloc_operator(opscr(1)%op)
        call dealloc_operator(opscr(2)%op)
        deallocate(opscr(1)%op)
        deallocate(opscr(2)%op)
        deallocate(opscr)
c      else
c        ! no external lines, so:
c        occ0 = 0
c        if(trim(name).ne.op_x_inter)then
c          call set_uop(op,name,.false.,
c     &         occ0,1,orb_info)
c        else
c          allocate(opscr(2))
c          allocate(opscr(1)%op)
c          allocate(opscr(2)%op)
c          call set_uop(opscr(1)%op,'scr1',.false.,
c     &         occ0,1,orb_info)
c
c          occ0(1,1:2) = 1
c          call set_uop(opscr(2)%op,'scr2',.false.,
c     &         occ0,1,orb_info)
c
c          call set_gen_intermediate(op,name,
c     &         opscr,2,orb_info)
c
c          call dealloc_operator(opscr(1)%op)
c          call dealloc_operator(opscr(2)%op)
c          deallocate(opscr(1)%op)
c          deallocate(opscr(2)%op)
c          deallocate(opscr)
c        endif
        do idx = 1, op%n_occ_cls
          op%formal_blk(idx) = .false.
        enddo

c      endif

      op%formal = .true.
      do idx = 1, op%n_occ_cls
        op%formal = op%formal.and.op%formal_blk(idx)
      end do

c      ! Extra terms required for CC-R12 intermediates.
c      ! Initially only for V.
c      if((trim(name).eq.op_v_inter.or.trim(name).eq.op_vbar_inter)
c     &     .and.cclog)then
c        if(trim(name).eq.op_v_inter)then
c          i=1
c          j=2
c        else
c          i=2
c          j=1
c        endif
c        allocate(opscr(2))
c        do idx = 1, 2
c          ! Shape of the operator.
c          occ0 = 0
c          call add_operator(op_dummy1,op_info)
c          idx_dum = idx_oplist2(op_dummy1,op_info)
c          dum1_pnt => op_info%op_arr(idx_dum)%op
c          occ0(2,i) = idx
c          call set_uop(dum1_pnt,op_dummy1,.false.,
c     &         occ0,1,orb_info)
c
cc          if(.not.r12fix)then
c            ! Spacer operator.
c            occ0 = 0
c            call add_operator(op_dummy2,op_info)
c            idx_dum = idx_oplist2(op_dummy2,op_info)
c            dum2_pnt => op_info%op_arr(idx_dum)%op
c            occ0(1,i) = 2
c            occ0(1,j) = 2-idx
c            call set_uop(dum2_pnt,op_dummy2,.false.,
c     &           occ0,1,orb_info)
c
c            opscr(1)%op => dum1_pnt
c            opscr(2)%op => dum2_pnt
c
c            ! Temporary intermediate.
c            call add_operator(op_dummy3,op_info)
c            idx_dum = idx_oplist2(op_dummy3,op_info)
c            dum3_pnt => op_info%op_arr(idx_dum)%op
c
c            call set_gen_intermediate(dum3_pnt,op_dummy3,
c     &           opscr,2,orb_info)
c
c            call join_operator(op,dum3_pnt,orb_info)
c
c            call del_operator(op_dummy3,op_info)
c            call del_operator(op_dummy2,op_info)
cc          else
cc            call join_operator(op,dum1_pnt,orb_info)
cc          endif
c          call del_operator(op_dummy1,op_info)
c        enddo
c        deallocate(opscr)
c      endif

      if (ntest.ge.100) then
        write(lulog,*) 'generated: '
        call print_op_occ(lulog,op)
        write(lulog,*) 'formal: ',op%formal_blk(1:op%n_occ_cls)
      end if

      return

      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public