MODULE module_sf_pxsfclay

 REAL , PARAMETER :: RICRIT = 0.25
 REAL , PARAMETER :: BETAH = 5.0
 REAL , PARAMETER :: BETAM = 5.0
 REAL , PARAMETER :: BM = 13.0
 REAL , PARAMETER :: BH = 15.7
 REAL , PARAMETER :: GAMAM = 19.3
 REAL , PARAMETER :: GAMAH = 11.6
 REAL , PARAMETER :: PR0 = 0.95
 REAL , PARAMETER :: CZO = 0.032
 REAL , PARAMETER :: OZO = 1.E-4
 REAL , PARAMETER :: VCONVC = 1.0

CONTAINS

   SUBROUTINE PXSFCLAY(U3D,V3D,T3D,TH3D,QV3D,P3D,dz8w, &
                     CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, &
                     ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
                     XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
                     U10,V10, &
                     GZ1OZ0,WSPD,BR,ISFFLX,DX, &
                     SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN, &
                     ids,ide, jds,jde, kds,kde, &
                     ims,ime, jms,jme, kms,kme, &
                     its,ite, jts,jte, kts,kte )

      IMPLICIT NONE
      INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte
      INTEGER, INTENT(IN ) :: ISFFLX
      REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
      REAL, INTENT(IN ) :: EP1,EP2,KARMAN
      REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
                INTENT(IN ) :: dz8w
      REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
                INTENT(IN ) :: QV3D, &
                                                              P3D, &
                                                              T3D, &
                                                             TH3D
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(IN ) :: MAVAIL, &
                                                             PBLH, &
                                                            XLAND, &
                                                              TSK
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(OUT ) :: U10, &
                                                              V10, &
                                                             QSFC
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: REGIME, &
                                                              HFX, &
                                                              QFX, &
                                                               LH, &
                                                          MOL,RMOL
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
                                                        PSIM,PSIH
      REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
                INTENT(IN ) :: U3D, &
                                                              V3D
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(IN ) :: PSFC
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: ZNT, &
                                                              ZOL, &
                                                              UST, &
                                                              CPM, &
                                                             CHS2, &
                                                             CQS2, &
                                                              CHS
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: FLHC,FLQC
      REAL, DIMENSION( ims:ime, jms:jme ) , &
                INTENT(INOUT) :: &
                                                              QGH
      REAL, INTENT(IN ) :: CP,G,ROVCP,R,XLV,DX
      REAL, DIMENSION( its:ite ) :: U1D, &
                                                              V1D, &
                                                             QV1D, &
                                                              P1D, &
                                                              T1D, &
                                                             TH1D
      REAL, DIMENSION( its:ite ) :: dz8w1d
      INTEGER :: I,J
      DO J=jts,jte
        DO i=its,ite
          dz8w1d(i) =dz8w(i,1,j)
          U1D(i) =U3D(i,1,j)
          V1D(i) =V3D(i,1,j)
          QV1D(i) =QV3D(i,1,j)
          P1D(i) =P3D(i,1,j)
          T1D(i) =T3D(i,1,j)
          TH1D(i) =TH3D(i,1,j)
        ENDDO
        CALL PXSFCLAY1D(J,U1D,V1D,T1D,TH1D,QV1D,P1D,dz8w1d, &
                CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j), &
                CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j), &
                ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j), &
                MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j), &
                XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j), &
                U10(ims,j),V10(ims,j), &
                FLHC(ims,j),FLQC(ims,j),QGH(ims,j), &
                QSFC(ims,j),LH(ims,j), &
                GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX, &
                SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN, &
                ids,ide, jds,jde, kds,kde, &
                ims,ime, jms,jme, kms,kme, &
                its,ite, jts,jte, kts,kte )
      ENDDO
   END SUBROUTINE PXSFCLAY
   SUBROUTINE PXSFCLAY1D(J,US,VS,T1D,THETA1,QV1D,P1D,dz8w1d, &
                     CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
                     ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
                     XLAND,HFX,QFX,TG, &
                     U10,V10,FLHC,FLQC,QGH, &
                     QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX, &
                     SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN, &
                     ids,ide, jds,jde, kds,kde, &
                     ims,ime, jms,jme, kms,kme, &
                     its,ite, jts,jte, kts,kte )
      IMPLICIT NONE
      REAL, PARAMETER :: XKA=2.4E-5
      REAL, PARAMETER :: PRT=1.
      INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte, &
                                        J
      INTEGER, INTENT(IN ) :: ISFFLX
      REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
      REAL, INTENT(IN ) :: EP1,EP2,KARMAN
      REAL, DIMENSION( ims:ime ) , &
                INTENT(IN ) :: MAVAIL, &
                                                             PBLH, &
                                                            XLAND, &
                                                              TG
      REAL, DIMENSION( ims:ime ) , &
                INTENT(IN ) :: PSFCPA
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: REGIME, &
                                                              HFX, &
                                                              QFX, &
                                                         MOL,RMOL
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
                                                        PSIM,PSIH
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: ZNT, &
                                                              ZOL, &
                                                              UST, &
                                                              CPM, &
                                                             CHS2, &
                                                             CQS2, &
                                                              CHS
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: FLHC,FLQC
      REAL, DIMENSION( ims:ime ) , &
                INTENT(INOUT) :: &
                                                              QGH
      REAL, DIMENSION( ims:ime ) , &
                INTENT(OUT) :: U10,V10, &
                                                          QSFC,LH
      REAL, INTENT(IN ) :: CP,G,ROVCP,XLV,DX,R
      REAL, DIMENSION( its:ite ), INTENT(IN ) :: dz8w1d
      REAL, DIMENSION( its:ite ), INTENT(IN ) :: US, &
                                                               VS, &
                                                             QV1D, &
                                                              P1D, &
                                                              T1D, &
                                                           THETA1
      REAL, DIMENSION( its:ite ) :: ZA, &
                                                              TH0, &
                                                           THETAG, &
                                                               WS, &
                                                            RICUT, &
                                                             USTM, &
                                                               RA, &
                                                          THETAV1, &
                                                         MOLENGTH
      REAL, DIMENSION( its:ite ) :: &
                                                      RHOX,GOVRTH
      REAL, DIMENSION( its:ite ) :: PSFC
      INTEGER :: KL
      INTEGER :: N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
      REAL :: PL,THCON,TVCON,E1
      REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
      REAL :: DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2
      REAL :: FLUXC,VSGD
      REAL :: XMOL,ZOBOL,Z10OL,ZNTOL,YNT,YOB,X1,X2
      REAL :: G2OZ0,G10OZ0,RA2,ZOLL
      REAL :: TV0,CPOT,RICRITI,AM,AH,SQLNZZ0,RBH,RBW,TSTV
      REAL :: PSIH2, PSIM2, PSIH10, PSIM10, CQS
      DO i = its,ite
        PSFC(I) = PSFCPA(I)/1000.
        TVCON = 1.0 + EP1 * QV1D(I)
        THETAV1(I)= THETA1(I) * TVCON
        RHOX(I) = PSFCPA(I)/(R*T1D(I)*TVCON)
      ENDDO
      DO I=its,ite
        E1=SVP1*EXP( SVP2*(TG(I)-SVPT0)/(TG(I)-SVP3) )
        QSFC(I)=EP2*E1/(PSFC(I)-E1)
        E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))
        PL = P1D(I)/1000.
        QGH(I)=EP2*E1/(PL-E1)
        CPM(I)=CP*(1.+0.8*QV1D(I))
      ENDDO
      DO I = its, ite
        TV0 = TG(I) * (1.0 + EP1 * QSFC(I)*MAVAIL(I))
        CPOT = (100./PSFC(I))**ROVCP
        TH0(I) = TV0 * (100./PSFC(I))**ROVCP
        THETAG(I) = CPOT * TG(I)
      ENDDO
      DO I = its,ite
        ZA(I) = 0.5 * DZ8W1D(I)
        WS(I) = SQRT(US(I) * US(I) + VS(I) * VS(I))
      ENDDO
        RICRITI = 1.0 / RICRIT
      DO i = its,ite
        GZ1OZ0(I) = ALOG(ZA(I) / ZNT(I))
        DTHVDZ = THETAV1(I) - TH0(I)
        fluxc = max(hfx(i)/rhox(i)/cp &
              + ep1*TH0(I)*qfx(i)/rhox(i),0.)
        VCONV = vconvc*(g/tg(i)*pblh(i)*fluxc)**.33
        VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
        WSPD(I)=SQRT(WS(I)*WS(I)+VCONV*VCONV+vsgd*vsgd)
        WSPD(I) = AMAX1(WSPD(I),0.1)
        GOVRTH(I) = G / THETA1(I)
        BR(I) = GOVRTH(I) * ZA(I) * DTHVDZ / (WSPD(I) * WSPD(I))
        RICUT(I) = 1.0 / (RICRITI + GZ1OZ0(I))
      ENDDO
      DO I = its,ite
        ZOLL = 0.0
        IF (BR(I) .GE. RICUT(I)) THEN
            REGIME(I) = 1.0
            ZOLL = BR(I) * GZ1OZ0(I) / (1.0 - RICRITI * RICUT(I))
            PSIM(I) = 1.0 - BETAM - ZOLL
            PSIH(I) = 1.0 - BETAH - ZOLL
        ELSE IF (BR(I) .GE. 0.0) THEN
            REGIME(I) = 2.0
            ZOLL = BR(I) * GZ1OZ0(I) / (1.0 - RICRITI * BR(I))
            PSIM(I) = -BETAM * ZOLL
            PSIH(I) = -BETAH * ZOLL
        ELSE
            REGIME(I) = 3.0
            AM = 0.031 + 0.276 * ALOG(GZ1OZ0(I))
            AH = 0.04 + 0.355 * ALOG(GZ1OZ0(I))
            SQLNZZ0 = SQRT(GZ1OZ0(I))
            PSIM(I) = AM * ALOG(1.0 - BM * SQLNZZ0 * BR(I))
            PSIH(I) = AH * ALOG(1.0 - BH * SQLNZZ0 * BR(I))
        ENDIF
      ENDDO
      DO I = its,ite
        DTG = THETA1(I) - THETAG(I)
        PSIX = GZ1OZ0(I) - PSIM(I)
        UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX
        USTM(I) = UST(I)
        IF ((XLAND(I)-1.5) .GE. 0.0) THEN
          ZNT(I) = CZO * USTM(I) * USTM(I) / G + OZO
          GZ1OZ0(I) = ALOG(ZA(I) / ZNT(I))
          PSIX = GZ1OZ0(I) - PSIM(I)
          UST(I) = KARMAN * WSPD(I) / PSIX
          USTM(I) = UST(I)
        ENDIF
        RA(I) = PR0 * (GZ1OZ0(I) - PSIH(I)) / (KARMAN * UST(I))
        RBH = 5.0 / UST(I)
        RBW = 4.47/UST(I)
        CHS(I) = 1./(RA(I) + RBH)
        CQS = 1./(RA(I) + RBW)
        MOL(I) = DTG * CHS(I) / UST(I)
        TSTV = (THETAV1(I) - TH0(I)) * CHS(I) / UST(I)
        IF (ABS(TSTV) .LT. 1.E-5) TSTV = 1.E-5
        MOLENGTH(I) = THETAV1(I) * UST(I) * UST(I) / (KARMAN * &
                        G * TSTV)
        XMOL = MOLENGTH(I)
        IF(MOLENGTH(I).GT.0.0) XMOL = AMAX1(MOLENGTH(I),2.0)
        RMOL(I) = 1/XMOL
        ZOL(I) = ZA(I)*RMOL(I)
        ZOBOL = 1.5*RMOL(I)
        Z10OL = 10.0*RMOL(I)
        ZNTOL = ZNT(I)*RMOL(I)
        IF(XMOL.LT.0.0) THEN
          YNT = ( 1.0 - GAMAH * ZNTOL )**0.5
          YOB = ( 1.0 - GAMAH * ZOBOL )**0.5
          PSIH2 = 2. * ALOG((YOB+1.0)/(YNT+1.0))
          x1 = (1.0 - gamam * z10ol)**0.25
          x2 = (1.0 - gamam * zntol)**0.25
          psim10 = 2.0 * ALOG( (1.0+x1) / (1.0+x2) ) + &
                         ALOG( (1.0+x1*x1) / (1.0+x2*x2)) - &
                         2.0 * ATAN(x1) + 2.0 * ATAN(x2)
        ELSE
          IF((ZOBOL-ZNTOL).LE.1.0) THEN
            PSIH2 = -BETAH*(ZOBOL-ZNTOL)
          ELSE
            PSIH2 = 1.-BETAH-(ZOBOL-ZNTOL)
          ENDIF
          IF((Z10OL-ZNTOL).LE.1.0) THEN
            PSIM10 = -BETAM*(Z10OL-ZNTOL)
          ELSE
            PSIM10 = 1.-BETAM-(Z10OL-ZNTOL)
          ENDIF
        ENDIF
        G2OZ0 = ALOG(1.5 / ZNT(I))
        G10OZ0 = ALOG(10.0 / ZNT(I))
        RA2 = PR0 * (G2OZ0 - PSIH2) / (KARMAN * UST(I))
        CHS2(I) = 1.0/(RA2 + RBH)
        CQS2(I) = 1.0/(RA2 + RBW)
        U10(I) = US(I)*(G10OZ0-PSIM10)/PSIX
        V10(I) = VS(I)*(G10OZ0-PSIM10)/PSIX
        FLHC(i)=CPM(I)*RHOX(I)*CHS(I)
        FLQC(i)=RHOX(I)*CQS*MAVAIL(I)
        QFX(I)=FLQC(I)*(QSFC(I)-QV1D(I))
        QFX(I)=AMAX1(QFX(I),0.)
        LH(I)=XLV*QFX(I)
        IF(XLAND(I)-1.5.GT.0.)THEN
          HFX(I)=-FLHC(I)*DTG
        ELSEIF(XLAND(I)-1.5.LT.0.)THEN
          HFX(I)=-FLHC(I)*DTG
          HFX(I)=AMAX1(HFX(I),-250.)
        ENDIF
      ENDDO
   END SUBROUTINE PXSFCLAY1D
   SUBROUTINE pxsfclayinit( allowed_to_read )
   LOGICAL , INTENT(IN) :: allowed_to_read
   INTEGER :: N
   REAL :: ZOLN,X,Y
   END SUBROUTINE pxsfclayinit
END MODULE module_sf_pxsfclay


text:
subroutine spline(n, x, y, b, c, d)
implicit none
integer n
real(8) x(n), y(n), b(n), c(n), d(n)
integer nm1, ib, i
real(8) t

nm1 = n-1
if(n<2) return
if(n<3) then
  b(1) = (y(2)-y(1))/(x(2)-x(1))
  c(1) = 0.
  d(1) = 0.
  b(2) = b(1)
  c(2) = 0.
  d(2) = 0.
  return
end if

d(1) = x(2) - x(1)
c(2) = (y(2) - y(1))/d(1)
do i = 2, nm1
  d(i) = x(i+1) - x(i)
  b(i) = 2.*(d(i-1) + d(i))
  c(i+1) = (y(i+1) - y(i))/d(i)
  c(i) = c(i+1) - c(i)
end do

b(1) = -d(1)
b(n) = -d(n-1)
c(1) = 0.
c(n) = 0.
if(n>3) then
  c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1))
  c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3))
  c(1) = c(1)*d(1)**2/(x(4)-x(1))
  c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3))
end if

do i = 2, n
  t = d(i-1)/b(i-1)
  b(i) = b(i) - t*d(i-1)
  c(i) = c(i) - t*c(i-1)
end do

c(n) = c(n)/b(n)
do ib = 1, nm1
  i = n-ib
  c(i) = (c(i) - d(i)*c(i+1))/b(i)
end do

b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2.*c(n))
do i = 1, nm1
  b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.*c(i))
  d(i) = (c(i+1) - c(i))/d(i)
  c(i) = 3.*c(i)
end do
c(n) = 3.*c(n)
d(n) = d(n-1)
return

end subroutine spline

real(8) function seval(n, u, x, y, b, c, d)
implicit none
integer n
real(8) u, x(n), y(n), b(n), c(n), d(n)
integer i, j, k
real(8) dx

i = 1
if(u>x(2)) then
  j = n+1
  do while(j>i+1)
    k = (i+j)/2
    if (u<x(k)) j = k
    if (u>=x(k)) i = k
  end do
end if

dx = u - x(i)
seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
return

end function seval

max_stars_repo_name: sdhzhs/ransfoil

max_issues_repo_name: sdhzhs/ransfoil

max_forks_repo_name: sdhzhs/ransfoil


text: array A = [3] string size 10 default "a"

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
program preprocessing
implicit none
integer, parameter :: n = 5
integer :: i

do i = 1, n
    print '(A, I0)', 'hello ', i
end do

#ifdef SAY_BYE
do i = 1, n
    print '(A, I0)', 'bye ', i
end do
#endif

end program preprocessing

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


text:
program main
    use MathClass
    use ArrayClass
    use LinearSolverClass
    implicit none

    integer(int32) :: i, num_data,sampling_rate
    real(real64),allocatable :: coeff(:),fit_coeff(:)
    real(real64),allocatable :: true_answer(:)
    real(real64),allocatable :: training_data(:,:)

    real(real64) :: err
    type(IO_) :: f
    type(Random_) :: random

    ! create training data
    true_answer = [0.5d0,100.0d0,0.10d0, 1.0d0, 200.0d0, -0.20d0]
    call f%open("step_response.txt","w")
    num_data = 50000
    sampling_rate = 100 !Hz

    do i=1,num_data
        write(f%fh,*) dble(i-1)/dble(sampling_rate) ,&
            offset_function(t=dble(i-1)/dble(sampling_rate) ,params=true_answer)&
            +random%gauss(mu=0.0d0,sigma=0.010d0)
    enddo

    call f%close()
    call f%plot(option="with lines")

    ! read training data
    training_data = zeros(num_data,2)
    call f%open("step_response.txt")
    do i=1,size(training_data,1)
        read(f%fh,*) training_data(i,1:2)
    enddo
    call f%close()

    ! initial guess
    !(1)coeff = [0.3d0,50.0d0,0.00d0, 1.0d0, 200.0d0, 0.00d0]
    coeff =          [0.10d0,  100.0d0,  0.0d0]
    coeff = coeff // [0.10d0,  200.0d0,  0.0d0]
    coeff = coeff // [0.10d0,  300.0d0,  0.0d0]
    coeff = coeff // [0.10d0,  400.0d0,  0.0d0]
    coeff = coeff // [0.10d0,  500.0d0,  0.0d0]
    coeff = coeff // [0.10d0,  600.0d0,  0.0d0]
    coeff = coeff // [0.10d0,  700.0d0,  0.0d0]

    ! fitting by Robbin-Monro algorithm (SGD)
    do i_i=1,10000

        fit_coeff = fit(f=offset_function,training_data=training_data, params=coeff ,&
            eta=0.00001d0,max_itr=30, error=err,use_ratio=0.0000030d0,logfile="SGD.log")
        coeff = fit_coeff
        print *, fit_coeff(:)

        call f%open("SGD_fit.txt","w")
        do i=0,num_data
            call f%write(dble(i)/100.0d0, offset_function(t=dble(i)/100.0d0,params=fit_coeff ) )
        enddo
        call f%close()

    enddo
    ! show result
    call f%open("training_data.txt","w")
    do i=1,size(training_data,1)
        call f%write(training_data(i,1), training_data(i,2) )
    enddo
    call f%close()

    print *, "--- check ---"

    print *, fit_coeff(:)
    print *, true_answer(:)

contains

real(real64) function offset_function(t,params)
    real(real64),intent(in) :: t
    real(real64),intent(in) :: params(:)
    integer(int32) :: n
    integer(int32) :: ii

    n = (size(params) - mod(size(params),3))/3
    offset_function = 0.0d0

    do ii=1,n
        offset_function = offset_function + sigmoid(x=t,params=params( (ii-1)*3+1:(ii-1)*3+3 ) )
    enddo

end function

end program main

max_stars_repo_name: kazulagi/faba

max_issues_repo_name: kazulagi/faba

max_forks_repo_name: kazulagi/faba


text:
!==============================================================================!
  subroutine User_Mod_End_Of_Time_Step(grid, n, time)
!------------------------------------------------------------------------------!
!   This function is called at the beginning of time step.                     !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Grid_Mod
  use Flow_Mod
  use Const_Mod, only: PI
  use Comm_Mod,  only: Comm_Mod_Global_Max_Real, this_proc
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
  integer         :: n     ! time step
  real            :: time  ! physical time
!----------------------------------[Locals]------------------------------------!
  integer :: c, e, dir
  real    :: lo, xo(4), yo(4), zo, ro, xc, yc, zc, vc, wc, sig_x, sig_yz,  &
             rmin, rmax, sg, lx, ly, lz, vmax
!==============================================================================!

  ! If not time for disturbing the velocity field, return
  if(mod(n-33, 1500) .ne. 0) return

  ! If too late to disturb, get out too
  if(n > 3000) return

  ! Print a message
  if(this_proc < 2) then
    print *, '# Superimposing random eddies on top of velocity field!'
  end if

  ! Minimum and maximum size of eddies
  rmin = 0.2
  rmax = 0.4

  ! Size of the computational domain
  lx = 2.0 * PI
  ly = PI
  lz = 2.0

  ! Browse through all eddies
  do e = 1, 80

    ! Random direction of the vortex
    call random_number(sg);
    if(sg < 0.5) then
      sg = -1.0
    else
      sg = +1.0
    end if

    ! Determine random position of a vortex
    call random_number(ro);     ro    = rmin + (rmax-rmin)*ro  ! rmin -> rmax
    call random_number(xo(1));  xo(1) = xo(1) * lx
    call random_number(yo(1));  yo(1) = yo(1) * ly
    call random_number(zo);     zo = ro + (lz - 2.0*ro) * zo

    ! Handle periodicity; that is: copy eddie in periodic directions
    xo(2:4) = xo(1)
    yo(2:4) = yo(1)
    if(xo(1) > lx/2.0) xo(3) = xo(1) - lx
    if(xo(1) < lx/2.0) xo(3) = xo(1) + lx
    if(yo(1) > ly/2.0) yo(2) = yo(1) - ly
    if(yo(1) < ly/2.0) yo(2) = yo(1) + ly
    xo(4) = xo(3)
    yo(4) = yo(2)

    ! Length of the eddy is three times
    lo = ro * 3.0

    sig_yz = ro / 2.0
    sig_x  = lo / 2.0

    ! Superimpose eddies on the velocity field
    do dir = 1, 4
      do c = 1, grid % n_cells
        xc = grid % xc(c)
        yc = grid % yc(c)
        zc = grid % zc(c)
        wc = sg * sin( (yc-yo(dir))/ro * PI )
        vc = sg * sin( (zc-zo)/ro * PI )

        vc = vc / (sig_yz*sqrt(PI+PI))*exp(-0.5*((yc-yo(dir))/sig_yz)**2)
        vc = vc / (sig_yz*sqrt(PI+PI))*exp(-0.5*((zc-zo)     /sig_yz)**2)

        wc = wc / (sig_yz*sqrt(PI+PI))*exp(-0.5*((yc-yo(dir))/sig_yz)**2)
        wc = wc / (sig_yz*sqrt(PI+PI))*exp(-0.5*((zc-zo)     /sig_yz)**2)

        vc = vc / (sig_x *sqrt(PI+PI))*exp(-0.5*((xc-xo(dir))/sig_x)**2)
        wc = wc / (sig_x *sqrt(PI+PI))*exp(-0.5*((xc-xo(dir))/sig_x)**2)

        v % n(c) = v % n(c) + vc
        w % n(c) = w % n(c) + wc
      end do
    end do
  end do

  vmax = 0
  do c = 1, grid % n_cells
    vmax = max(vmax, abs(v % n(c)))
    vmax = max(vmax, abs(w % n(c)))
  end do
  call Comm_Mod_Global_Max_Real(vmax)
  do c = 1, grid % n_cells
    v % n(c) = v % n(c) / vmax / 10.0
    w % n(c) = w % n(c) / vmax / 10.0
  end do

  end subroutine

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


text:
*DK ung2avs
      subroutine ung2avs(imsgin,xmsgin,cmsgin,msgtype,nwds,ics)
C
C#######################################################################
C
C      PURPOSE -
C
C      This routine is used to convert files in UNGenerate format
C      to AVS format. It allows the user to input a constant value for
C      the third coordinate. The default is 0.0. The xyz can be
C      reordered in any combination when written to the AVS file by
C      defining the two columns of coordinate date in the ung2avs
C      command to be xy, xz, yx, zx, yz, or zy. The default is xy.
C
C          FORMAT: UNG2AVS/avs_file_out/ung_file_in/[constant]/ &
C                                      [xy, xz, yx, zx, yz, zy]
C
C                 default is:
C                 UNG2AVS/avs_file_out/ung_file_in
C                     converts ung_file_in to avs_file_out
C                     with a constant default z value of 0.0
C                     in xy order
C                       or
C                 UNG2AVS/avs_file_out/ung_file_in/7.0
C                     converts ung_file_in to avs_file_out
C                     with a constant default z value of 7.0
C                     in xy order
C                       or
C                 UNG2AVS/avs_file_out/ung_file_in/15.0/xz
C                     converts ung_file_in to avs_file_out
C                     with a constant y value of 15.0
C                     in xz order
C
C          UNG FILE FORMAT:
C              The UNG format contains a file of xy values, in groups of
C              sequentially numbered sets, of connected line segments,
C              each ending with END. The format for connected lines that
C              form polygons or that do not form polygons is very similar.
C              To make sets of connected lines that are polygons the first
C              xy pair for each set of lines must be on the same line as
C              the set number. Closed polygons will be formed by connecting
C              the last point of a point set to the first point of the
C              point set. The set number is on a line by itself for
C              non-polygons. Spacing on a line does not matter.
C
C              FORMAT for connected lines that do not form polygons:
C
C                     1
C                       x   y
C                       x   y
C                       x   y
C                         ...
C                     END
C                     2
C                       x   Y
C                       x   y
C                       x   y
C                         ...
C                     END
C                         ...
C                     n
C                       x   y
C                       x   y
C                         ...
C                     END
C                     END
C
C              FORMAT for polygons:
C
C                     1 x   y
C                       x   y
C                       x   y
C                         ...
C                     END
C                     2 x   y
C                       x   y
C                       x   y
C                         ...
C                     END
C                         ...
C                     n x   y
C                       x   y
C                         ...
C                     END
C                     END
C
C              SAMPLE of connected, non-polygon, lines:
C
C                  UNG file                          AVS file
C                1                           10  7 0 0 0
C                  10.00  15.00               1   10.0000 15.0000  0.0000
C                  10.50  15.50               2   10.5000 15.5000  0.0000
C                  11.00  16.00               3   11.0000 16.0000  0.0000
C                END                          4   20.0000 25.0000  0.0000
C                2                            5   20.5000 25.5000  0.0000
C                  20.00  25.00               6   21.0000 26.0000  0.0000
C                  20.50  25.50               7   22.0000 27.0000  0.0000
C                  21.00  26.00               8   30.0000 35.0000  0.0000
C                  22.00  27.00               9   30.5000 35.5000  0.0000
C                END                         10   31.0000 36.0000  0.0000
C                3                           1  1 line    1    2
C                  30.00  35.00              2  1 line    2    3
C                  30.50  35.50              3  2 line    4    5
C                  31.00  36.00              4  2 line    5    6
C                END                         5  2 line    6    7
C                END                         6  3 line    8    9
C                                            7  3 line    9   10
C
C              SAMPLE of polygons:
C
C                  UNG file                          AVS file
C                1 10.00  15.00              10  10  0 0 0
C                  10.50  15.50               1   10.0000 15.0000  0.0000
C                  11.00  16.00               2   10.5000 15.5000  0.0000
C                END                          3   11.0000 16.0000  0.0000
C                2 20.00  25.00               4   20.0000 25.0000  0.0000
C                  20.50  25.50               5   20.5000 25.5000  0.0000
C                  21.00  26.00               6   21.0000 26.0000  0.0000
C                  22.00  27.00               7   22.0000 27.0000  0.0000
C                END                          8   30.0000 35.0000  0.0000
C                3 30.00  35.00               9   30.5000 35.5000  0.0000
C                  30.50  35.50              10   31.0000 36.0000  0.0000
C                  31.00  36.00              1  1 line    1    2
C                END                         2  1 line    2    3
C                END                         3  1 line    3    1
C                                            4  2 line    4    5
C                                            5  2 line    5    6
C                                            6  2 line    6    7
C                                            7  2 line    7    4
C                                            8  3 line    8    9
C                                            9  3 line    9   10
C                                           10  3 line   10    8
C
C      INPUT ARGUMENTS -
C
C         NONE
C
C
C      OUTPUT ARGUMENTS -
C
C         NONE
C
C
C      CHANGE HISTORY -
C
C         $Log: ung2avs.f,v $
C         Revision 2.00  2007/11/09 20:04:05  spchu
C         Import to CVS
C
CPVCS
CPVCS       Rev 1.6   18 May 2006 15:04:56   gable
CPVCS    Modified output format so code will work with more than 9999 elements.
CPVCS
CPVCS       Rev 1.5   26 Feb 2004 10:58:04   dcg
CPVCS    use explicit size for temporary arrays - linux will
CPVCS    not compile with nwds as size
CPVCS
CPVCS       Rev 1.4   26 Nov 2003 11:41:34   gable
CPVCS    Changed coordinate arrays to double precision arrays.
CPVCS    Modified output format statement to more significant figures
CPVCS    and changed some integer output formats so they would work
CPVCS    with bigger integers.
CPVCS
CPVCS       Rev 1.3   Mon Dec 13 12:49:00 1999   llt
CPVCS    modified documentation
CPVCS
CPVCS       Rev 1.2   Wed Dec 08 13:39:40 1999   llt
CPVCS    added option to allow for all combinations of the xyz ordering
CPVCS
CPVCS       Rev 1.1   Fri Jan 22 14:43:56 1999   llt
CPVCS    modified header
CPVCS
CPVCS       Rev 1.0   Fri Jan 22 14:28:16 1999   llt
CPVCS    Initial revision.
C
C#######################################################################
C
      implicit none
c
c      include "local_element.h"
C
C#######################################################################
C
c parser variables for subroutine call
      integer nwds, imsgin(nwds), msgtype(nwds)
      real*8 xmsgin(nwds)
      character*(*) cmsgin(nwds)

      integer ics, ierror, icscode, iwrite
      character*40 isubname
c
      integer iunit_in, iunit_out
      logical trueend, polytrue, firstmaterial
      integer nnodes, nmat, i, mat_num, j, i_save, length
      pointer (ipmat_nums, mat_nums(1000000))
      integer mat_nums
      real*8 z_value, x_save
      pointer (ipx_value, x_value(1000000))
      pointer (ipy_value, y_value(1000000))
      real*8 x_value, y_value
      pointer (ipline_end, line_end(10000000))
      logical line_end
      character*2 order

c local parser variables
      character*250 input_msg
      integer lenparse
      integer msg(10)
      real*8 xmsg(10)
      integer imsg(10)
      character*32 cmsg(10), file_avs, file_ung
C
C
C#######################################################################
C
      character*132 logmess
C
C#######################################################################
C
      isubname='ung2avs'
      iunit_in = -1
      iunit_out = -1
      ierror = 0
      iwrite = 1
C
      if(nwds .eq. 3 .or. nwds .eq. 4 .or. nwds .eq. 5) then

        file_avs =  cmsgin(2)
        call hassign (iunit_out, cmsgin(2), ierror)
        if (iunit_out.lt.0 .or. ierror.lt.0) then
          iwrite = 0
          ierror = 1
          call x3d_error(isubname,'out hassign bad file unit')
          goto 9999
        endif

        file_ung =  cmsgin(3)
        call hassign (iunit_in, cmsgin(3), ierror)
        if (iunit_in.lt.0 .or. ierror.lt.0) then
          iwrite = 0
          ierror = 1
          call x3d_error(isubname,'in  hassign bad file unit')
          goto 9999
        endif

        if (nwds .eq. 4 .or. nwds .eq. 5) then
          z_value = xmsgin(4)
        else
          z_value = 0.000
        end if
        order = 'xy'
        if (nwds .eq. 5) then
          if (cmsgin(5) .eq. 'xy') then
            order = 'xy'
          elseif (cmsgin(5) .eq. 'xz') then
            order = 'xz'
          elseif (cmsgin(5) .eq. 'yx') then
            order = 'yx'
          elseif (cmsgin(5) .eq. 'zx') then
            order = 'zx'
          elseif (cmsgin(5) .eq. 'yz') then
            order = 'yz'
          elseif (cmsgin(5) .eq. 'zy') then
            order = 'zy'
          else
            order = 'xy'
          end if
        end if

      else
         write(logmess,'(a)') "ERROR: command not recognized."
         call writloga('default',0,logmess,0,ics)
         write(logmess,3000)
 3000    format(' UNG2AVS / file_out / file_in / [z_value] ')
         call writloga('default',0,logmess,0,ics)
         goto 9999
      endif
C
c initialize values
c
      trueend = .false.
      nnodes = 0
      nmat = 0
      firstmaterial = .true.
      polytrue = .false.
      nwds = 0

c     read file to calculate maximum number of nodes -- for memory management
c     change loop to avoid length restriction and check for empty file
      length = 0
      do
        read (iunit_in,*, END=90)
        length = length + 1
      end do
90    rewind (iunit_in)

      if (length .le. 0) then
         write(logmess,'(a,a)') "ERROR: empty file: ",file_ung
         call writloga('default',0,logmess,0,ics)
         ierror = ierror + 1
         iwrite = 0
         goto 9999
      endif

c allocate memory

      call mmgetblk ('mat_nums',isubname,ipmat_nums,length,2,icscode)
      call mmgetblk ('x_value',isubname,ipx_value,length,2,icscode)
      call mmgetblk ('y_value',isubname,ipy_value,length,2,icscode)
      call mmgetblk ('line_end',isubname,ipline_end,length,2,icscode)

c read and parse string
c      do i=1, 1000000
       do i=1, length
         read(iunit_in,'(a)', end=100) input_msg
         lenparse = len(input_msg)

         call parse_string2(lenparse, input_msg,
     .                     imsg,msg,xmsg,cmsg,nwds)

c check to see if at a new material
         if ((nwds .eq. 1 .or. nwds .eq. 3) .and.
     .        (cmsg(1)(1:3) .ne. 'END')) then
           trueend = .false.
           mat_num = imsg(1)
           nmat = nmat + 1
           if (nwds .eq. 3) then
             polytrue = .true.
             nnodes = nnodes + 1
             mat_nums(nnodes) = mat_num
             x_value(nnodes) = xmsg(2)
             y_value(nnodes) = xmsg(3)
             line_end(nnodes) = .false.
           endif
c check to see if at END
         else if ((nwds .eq. 1) .and.
     .            (cmsg(1)(1:3) .eq. 'END') .and.
     .            (.not. trueend)) then
           trueend = .true.
           line_end(nnodes) = .true.
c check to see if at END END
         else if ((nwds .eq. 1) .and. (cmsg(1)(1:3) .eq. 'END') .and.
     .            trueend) then
           goto 100
c calculate number of nodes and write x,y
         elseif (nwds .eq. 2) then
            if (polytrue) then
              if (firstmaterial) then
                x_save = xmsg(1)
                firstmaterial = .false.
              endif
            endif
            nnodes = nnodes + 1
            mat_nums(nnodes) = mat_num
            x_value(nnodes) = xmsg(1)
            y_value(nnodes) = xmsg(2)
            line_end(nnodes) = .false.
         end if
      end do

 100  continue
c

c WRITE AVS FILE
c
c   write header
      if (polytrue) then
        write (iunit_out, *) nnodes, nnodes, ' 0 0 0'
      else
        write (iunit_out, *) nnodes, ' ', nnodes-nmat, ' 0 0 0 '
      end if

c   write xyz values

        if (order .eq. 'xy') then
          do i=1, nnodes
            write(iunit_out, 110) i, x_value(i), y_value(i), z_value
          end do
        elseif (order .eq. 'xz') then
          do i=1, nnodes
            write(iunit_out, 110) i, x_value(i), z_value, y_value(i)
          end do
        elseif (order .eq. 'yx') then
          do i=1, nnodes
            write(iunit_out, 110) i, y_value(i), x_value(i), z_value
          end do
        elseif (order .eq. 'zx') then
          do i=1, nnodes
            write(iunit_out, 110) i, y_value(i), z_value, x_value(i)
          end do
        elseif (order .eq. 'yz') then
          do i=1, nnodes
            write(iunit_out, 110) i, z_value, x_value(i), y_value(i)
          end do
        elseif (order .eq. 'zy') then
          do i=1, nnodes
            write(iunit_out, 110) i, z_value, y_value(i), x_value(i)
          end do
        end if

 110  format(i10, 2x, 3(1pe20.12))

c   write connectivity
      j = 1
      i_save = 1
      do i=1,nnodes
        if (.not. line_end(i)) then
          write(iunit_out, 120) j, mat_nums(i),
     .                        i,  i+1
 120      format (i10, 2x, i6, ' line ', i10, 2x, i10)
           j = j + 1
        else
          if (polytrue) then
            write(iunit_out, 120) j, mat_nums(i),
     .                          i, i_save
            i_save = i+1
            j = j + 1
          end if
        end if
      end do

c deallocate memory

 9999 continue

      call mmrelblk ('mat_nums',isubname,ipmat_nums,icscode)
      call mmrelblk ('x_value',isubname,ipx_value,icscode)
      call mmrelblk ('y_value',isubname,ipy_value,icscode)
      call mmrelblk ('line_end',isubname,ipline_end, icscode)

      if (iunit_in.gt.0) then
        close (iunit_in)
        write(logmess,'(a,a)') "UNG2AVS file read: ",file_ung
        call writloga('default',0,logmess,0,ics)
        write(logmess,'(a,i15)')
     &  "UNG2AVS read number lines: ",length
        call writloga('default',0,logmess,0,ics)
      else
        ierror = ierror + 1
      endif

      if (ierror .ne. 0) then
        write(logmess,'(a,a)') "UNG2AVS file NOT read: ",file_ung
        call writloga('default',0,logmess,0,ics)
      endif

      if (iunit_out.gt.0) then
        close (iunit_out)
        write(logmess,'(a,a)') "UNG2AVS file written: ",file_avs
        call writloga('default',0,logmess,1,ics)
      else
        ierror = ierror + 1
      endif

      if (iwrite .ne. 1) then
        write(logmess,'(a,a)') "UNG2AVS file NOT written: ",file_avs
        call writloga('default',0,logmess,1,ics)
      endif

      return
      end

ext: f

max_forks_repo_name: daniellivingston/LaGriT


text:
subroutine do_nothing
  continue
end subroutine do_nothing

max_stars_repo_name: EraYaN/transpyle

max_issues_repo_name: EraYaN/transpyle

max_forks_repo_name: EraYaN/transpyle


text:
C    @(#)deleterpt.f	20.6 7/18/96
C****************************************************************
C
C   	File: deleterpt.f
C
C   	Purpose: Generates a report of all system deletions.
C
C   	Author: Walt Powell            Date: 13 November 1992
C   	Called by: p_report.f
C
C****************************************************************
C
      	integer function deleterpt (in_buffer, out_buffer, scrfil)
        character in_buffer *(*), out_buffer *(*)
        integer scrfil

      	include 'ipfinc/parametr.inc'

      	include 'ipfinc/blank.inc'
      	include 'ipfinc/changr.inc'
      	include 'ipfinc/lfiles.inc'
      	include 'ipfinc/delete.inc'

	character null * 1, linefeed * 1, header(4) * 80
        integer o2, apdoutbuf

        maxbuf_out = len( out_buffer ) - 400
        null = char(0)
        linefeed = char(10)
        deleterpt = 0

        write (header(1), 100) ndelete, cspare(30), dte

  100   format (' List of system deletions - Total ', i4, ' Case ', a,
     &     ' date ', a)
        header(2) = ' '

        out_buffer(1:1) = null
        o2 = index (out_buffer,null)
        do i = 1, 2
           if (scrfil .gt. 0) write (scrfil, '(a)') header(i)
           length = apdoutbuf(o2, header(i), out_buffer(o2:))
           o2 = o2 + length
        enddo
        ic = 1
        do while (ic .le. ndelete .and. o2 .le. maxbuf_out)
           if (delete(ic)(3:3) .ne. '*') then
              if (scrfil .gt. 0) write (scrfil, '(a)') delete(ic)(1:92)
              length = apdoutbuf(o2, delete(ic)(1:92), out_buffer(o2:))
              o2 = o2 + length
           endif
           ic = ic + 1
        enddo

        if (o2 .gt. maxbuf_out) then
           write (out_buffer(o2:o2+8), 820) linefeed, null
  820      format (a, '*[MORE]', a)
           o2 = o2 + 9
        endif

  900   continue
        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


module my_module
  interface my_interface
     module procedure my_sub
  end interface
contains
  subroutine my_sub
  end subroutine my_sub
end module my_module


text:
c
c-----------------------------------------------------------------------
c
c
c
      double precision function p18th(x)
      implicit double precision (a-h,o-z)
c
c     function specifying theta as a function of x along the lower wall
c     theta is the angle of the wall = angle of streamline
c
      if (x.le.0.1d0) slope = 0.d0
      if (x.gt.0.1d0 .and. x.lt.0.7d0) slope = -0.6d0*(x-0.1d0)
      if (x.ge.0.7d0) slope = -0.36d0
      p18th = datan(slope)
      return
      end


text:
module maths
use vars
! use Fast_Fourier
use singleton
!use nrutil, ONLY : assert_eq, imaxloc, outerprod, swap
implicit none
contains

!----------------------------------------------------------------
! This function calculates the 3-j symbol
! J_i and M_i have to be twice the actual value of J and M
!----------------------------------------------------------------
   function w3js(j1,j2,j3,m1,m2,m3)
   integer :: m1, m2, m3, j1, j2, j3
	integer :: ia, ib, ic, id, ie, im, ig, ih, z, zmin, zmax, jsum
	real(kind=8) :: w3js, cc, denom, cc1, cc2

!		w3js = w3js_regge(j1/2,j2/2,j3/2,m1/2,m2/2,m3/2)
!		return
      w3js = 0.d0
      if (m1+m2+m3 /= 0) return
      ia = j1 + j2
      if (j3 > ia) return
      ib = j1 - j2
      if (j3 < abs(ib)) return
		if (abs(m1) > j1) return
		if (abs(m2) > j2) return
		if (abs(m3) > j3) return

      jsum = j3 + ia
      ic = j1 - m1
      id = j2 - m2

		ie = j3 - j2 + m1
		im = j3 - j1 - m2
		zmin = max0(0,-ie,-im)
		ig = ia - j3
		ih = j2 + m2
		zmax = min0(ig,ih,ic)
		cc = 0.d0

		do z = zmin, zmax, 2
			denom = fact(z/2)*fact((ig-z)/2)*fact((ic-z)/2)*fact((ih-z)/2)*&
				fact((ie+z)/2)*fact((im+z)/2)
      	if (mod(z,4) /= 0) denom = -denom
			cc = cc + 1.d0/denom
		enddo

		cc1 = fact(ig/2)*fact((j3+ib)/2)*fact((j3-ib)/2)/fact((jsum+2)/2)
      cc2 = fact((j1+m1)/2)*fact(ic/2)*fact(ih/2)*fact(id/2)*fact((j3-m3)/2)*fact((j3+m3)/2)
      cc = cc * sqrt(1.d0*cc1*cc2)
		if (mod(ib-m3,4) /= 0) cc = -cc
		w3js = cc
		if (abs(w3js) < 1.d-8) w3js = 0.d0
1000 		return
   end function w3js

!----------------------------------------------------------------
! This function calculates the 6-j symbol
! J_i and M_i have to be twice the actual value of J and M
!----------------------------------------------------------------
	function w6js(j1,j2,j3,l1,l2,l3)
	integer :: j1,j2,j3,l1,l2,l3
	integer :: ia, ib, ic, id, ie, iif, ig, ih, sum1, sum2, sum3, sum4
	integer :: w, wmin, wmax, ii, ij, ik
   real(kind=8) :: w6js, omega, theta1, theta2, theta3, theta4, theta, denom

      w6js = 0.d0
		ia = j1 + j2
		if (ia < j3) return
		ib = j1 - j2
		if (abs(ib) > j3) return
		ic = j1 + l2
		if (ic < l3) return
		id = j1 - l2
		if (abs(id) > l3) return
		ie = l1 + j2
		if (ie < l3) return
		iif = l1 - j2
		if (abs(iif) > l3) return
		ig = l1 + l2
		if (ig < j3) return
		ih = l1 - l2
		if (abs(ih) > j3) return
      sum1=ia + j3
      sum2=ic + l3
      sum3=ie + l3
      sum4=ig + j3
		wmin = max0(sum1, sum2, sum3, sum4)
		ii = ia + ig
		ij = j2 + j3 + l2 + l3
		ik = j3 + j1 + l3 + l1
		wmax = min0(ii,ij,ik)
		omega = 0.d0
		do w = wmin, wmax, 2
	      denom = fact((w-sum1)/2)*fact((w-sum2)/2)*fact((w-sum3)/2)&
				*fact((w-sum4)/2)*fact((ii-w)/2)*fact((ij-w)/2)&
				*fact((ik-w)/2)
			if (mod(w,4) /= 0) denom = -denom
			omega = omega + fact(w/2+1) / denom
		enddo
      theta1 = fact((ia-j3)/2)*fact((j3+ib)/2)*fact((j3-ib)/2)/fact(sum1/2+1)
      theta2 = fact((ic-l3)/2)*fact((l3+id)/2)*fact((l3-id)/2)/fact(sum2/2+1)
      theta3 = fact((ie-l3)/2)*fact((l3+iif)/2)*fact((l3-iiF)/2)/fact(sum3/2+1)
      theta4 = fact((ig-j3)/2)*fact((j3+ih)/2)*fact((j3-ih)/2)/fact(sum4/2+1)
      theta = theta1 * theta2 * theta3 * theta4
		w6js = omega * sqrt(theta)
      if (abs(w6js) < 1.d-8) w6js = 0.d0
1000 		return
   end function w6js

!----------------------------------------------------------------
! This function calculates the 9-j symbol
! J_i and M_i have to be twice the actual value of J and M
!----------------------------------------------------------------
   function w9js(j1,j2,j3,j4,j5,j6,j7,j8,j9)
	integer :: j1,j2,j3,j4,j5,j6,j7,j8,j9
	integer :: i, kmin, kmax, k
	real(kind=8) :: x, s, x1, x2, x3, w9js

      kmin = abs(j1-j9)
		kmax = j1 + j9
		i = abs(j4-j8)
		if (i > kmin) kmin = i
		i = j4 + j8
		if (i < kmax) kmax = i
		i = abs(j2-j6)
		if (i > kmin) kmin = i
		i = j2 + j6
		if (i < kmax) kmax = i
		x = 0.d0
		do k = kmin, kmax, 2
			s = 1.d0
			if (mod(k,2) /= 0) s = -1.d0
			x1 = w6js(j1,j9,k,j8,j4,j7)
      	x2 = w6js(j2,j6,k,j4,j8,j5)
      	x3 = w6js(j1,j9,k,j6,j2,j3)
      	x = x + s*x1*x2*x3*dfloat(k+1)
		enddo
      w9js = x
      return
   end function w9js

!----------------------------------------------------------------
! Reorder the Regge square to make it similar to eq. (2.11) in Rasch & Yu (SIAM J. Sci. Comput. Vol 25, n 4, 1416)
!----------------------------------------------------------------
   function reorder_regge(regge, factor)
	integer :: reorder_regge(3,3), regge(3,3), temp(3,3), loc_max(2), col_max, row_max, loc_min(2), col_min, row_min, t
	integer :: i, j, tot, minus1, factor

		tot = sum(regge(1,:))
		minus1 = (-1)**tot
		factor = 1

		temp = regge

		loc_max = maxloc(temp)
		row_max = loc_max(2)
		col_max = loc_max(1)

		loc_min = minloc(temp)
		row_min = loc_min(2)
		col_min = loc_min(1)

! Put the maximum and minimum in the same row
		if (col_min == col_max) then
			temp = transpose(regge)

			row_max = loc_max(1)
			col_max = loc_max(2)

			row_min = loc_min(1)
			col_min = loc_min(2)

		endif

		reorder_regge = temp
		if (row_min /= 1) then
			reorder_regge(:,1) = temp(:,row_min)
			reorder_regge(:,row_min) = temp(:,1)
			factor = factor * minus1
		endif
		temp = reorder_regge

		loc_min = minloc(temp)
		row_min = loc_min(2)
		col_min = loc_min(1)

! Put smallest one into R_11
		reorder_regge = temp
		if (col_min /= 1) then
			reorder_regge(1,:) = temp(col_min,:)
			reorder_regge(col_min,:) = temp(1,:)
			factor = factor * minus1
			temp = reorder_regge
		endif

! Put largest one into R_21
		loc_max = maxloc(temp)
		row_max = loc_max(2)
		col_max = loc_max(1)

		if (col_max /= 2) then
			reorder_regge = temp
			reorder_regge(2,:) = temp(col_max,:)
			reorder_regge(col_max,:) = temp(2,:)
			factor = factor * minus1
		endif

   end function reorder_regge

!----------------------------------------------------------------
! Calculate the 3j symbol by indexing a 1D array with all the 3j symbols Rasch & Yu (SIAM J. Sci. Comput. Vol 25, n 4, 1416)
!----------------------------------------------------------------
	function w3js_regge(j1,j2,j3,m1,m2,m3)
	integer :: j1, j2, j3, m1, m2, m3, regge(3,3), regge_order(3,3), L, X, T, B, S, c1, c2, c3, c4, c5, c, factor
	real(kind=8) :: w3js_regge

		w3js_regge = 0.d0
		if (abs(m1) > j1 .or. abs(m2) > j2 .or. abs(m3) > j3 .or. (m1+m2+m3) /= 0) return

		regge(1,1) = -j1+j2+j3
		regge(2,1) = j1-j2+j3
		regge(3,1) = j1+j2-j3

		regge(1,2) = j1-m1
		regge(2,2) = j2-m2
		regge(3,2) = j3-m3

		regge(1,3) = j1+m1
		regge(2,3) = j2+m2
		regge(3,3) = j3+m3

		regge_order = reorder_regge(regge,factor)

		L = regge_order(2,1)
		X = regge_order(1,2)
		T = regge_order(3,3)
		B = regge_order(2,2)
		S = regge_order(1,1)

		c1 = L*(24+L*(50+L*(35+L*(10+L))))/120
		c2 = X*(6+X*(11+X*(6+X)))/24
		c3 = T*(2+T*(3+T))/6
		c4 = B*(B+1)/2
		c5 = S
		c = c1+c2+c3+c4+c5+1

		w3js_regge = factor*threej(c)
	end function w3js_regge

!----------------------------------------------------------------
! Calculate the list of 3j symbols and store them in a 1D array
!----------------------------------------------------------------
	subroutine init_regge(Lmax)
	integer :: S, B, T, X, L, Lmax, cmax, c1, c2, c3, c4, c5, c, k, j1, j2, j3, m1, m2, m3

		cmax = Lmax*(274+Lmax*(225+Lmax*(85+Lmax*(15+Lmax))))/120+1

		allocate(threej(0:cmax))

		k = 0
		do L = 0, Lmax
			do X = 0, L
				do T = 0, X
					do B = 0, T
						do S = 0, B
							j1 = (B+L-T+X)
							j2 = (B+S-T+X)
							j3 = (L+S)
							m1 = (B+L-T-X)
							m2 = (-B+S-T+X)
							m3 = (T+T-L-S)

							c1 = L*(24+L*(50+L*(35+L*(10+L))))/120
							c2 = X*(6+X*(11+X*(6+X)))/24
							c3 = T*(2+T*(3+T))/6
							c4 = B*(B+1)/2
							c5 = S
							c = c1+c2+c3+c4+c5+1

							threej(c) = w3js(j1,j2,j3,m1,m2,m3)

							k = k + 1
						enddo
					enddo
				enddo
			enddo
		enddo

		open(unit=12,file='regge_3j.dat',action='write',form='unformatted')
		write(12) threej
		close(12)

	end subroutine init_regge

!----------------------------------------------------------------
! This function calculates the rotation matrix of order k
! The results are returned in dr(iq,iqp) and di(iq,iqp) representing
! the real and imaginary part of the rotation matrix respectevely.
! The angles alpha, beta and gamma are given in radians
!----------------------------------------------------------------
   subroutine rotmatk(k,alpha,beta,gamma,dr,di)
	integer :: k
	real(kind=8) :: alpha, beta, gamma
   real(kind=8) :: dr(-1:1,-1:1), di(-1:1,-1:1)
	integer :: iq, iqp, it, it1, it2, it3, itmax, itmin, ie1, ie2
	real(kind=8) :: a1, a2, ca, sa, cg, sg, bo2, cb2, sb2, redm, s, den

		dr = 0.d0
		di = 0.d0
		do iq = -k, k
			do iqp = -k, k
				a1 = alpha * dfloat(iq)
				a2 = gamma * dfloat(iqp)
				ca = dcos(a1)
				sa = dsin(a1)
				cg = dcos(a2)
				sg = dsin(a2)
				bo2 = beta / 2.d0
				cb2 = dcos(bo2)
				sb2 = dsin(bo2)
				it1 = k + iq
				it2 = k - iqp
				itmax = it1
				if (it2 < itmax) itmax = it2
				itmin = 0
				it3 = iq - iqp
				if (it3 > 0) itmin = it3
				redm = 0.d0
				if (beta /= 0.d0) then
					do it = itmin, itmax
						s = 1.d0
						if (mod(it,2) /= 0) s = -1.d0
						ie1 = 2*(k-it)+iq-iqp
						ie2 = 2*it+iqp-iq
						den = fact(it1-it) * fact(it2-it) * fact(it) * fact(it-it3)
						redm = redm + s*(cb2**ie1)*(sb2**ie2) / den
					enddo
					redm = redm*sqrt(fact(k+iq)*fact(k-iq)*fact(k+iqp)*fact(k-iqp))
				else
					if (iq == iqp) redm = 1.d0
				endif
				dr(iq,iqp) = redm*(ca*cg-sa*sg)
				di(iq,iqp) = -redm*(ca*sg+sa*cg)
			enddo
		enddo

	end subroutine rotmatk

!----------------------------------------------------------------
! This function calculates the reduced rotation matrix of order 0,1,2,...kmax
! The results are returned in complex format dr(K,iq,iqp)
! The angle beta is given in radians
! The maximum value of kmax is 6
!----------------------------------------------------------------
   subroutine reduced_matrix(kmax,beta,reduced)
	integer :: kmax, K, M, N, t, min_t, max_t
	real(kind=8) :: beta, cos_beta, sin_beta, suma, f1, f2
	real(kind=8) :: reduced(0:kmax,-kmax:kmax,-kmax:kmax)

		cos_beta = cos(beta/2.d0)
      sin_beta = sin(beta/2.d0)

		do K = 0, kmax
			do M = -K, K
				do N = -K, K
					max_t = min(K+M,K-N)
        			min_t = max(0,M-N)
					suma = 0.d0
					f1 = sqrt(fact(K+M)*fact(K-M)*fact(K+N)*fact(K-N))
					do t = min_t, max_t
						f2 = fact(K+M-t)*fact(K-N-t)*fact(t)*fact(t+N-M)
						suma = suma + (-1.d0)**t*cos_beta**(2*K+M-N-2*t)*sin_beta**(2*t-M+N) / f2
					enddo
					reduced(k,M,N) = f1 * suma
				enddo
			enddo
		enddo

	end subroutine reduced_matrix

!----------------------------------------------------------------
! This function calculates the rotation matrix of order 0,1,2,...kmax
! The results are returned in complex format dr(K,iq,iqp)
! The angles alpha, beta and gamma are given in radians
! The maximum value of kmax is 6
!----------------------------------------------------------------
   subroutine rotmatsu(kmax,alpha,beta,gamma,dr)
   integer :: kmax, K, M, N
	real(kind=8) :: alpha, beta, gamma, ca, sa, cg, sg
	complex(kind=8) :: dr(0:kmax,-kmax:kmax,-kmax:kmax), ii
	real(kind=8) :: reduced(0:kmax,-kmax:kmax,-kmax:kmax)

		ii = cmplx(0.d0,1.d0)
		call reduced_matrix(kmax,beta,reduced)

		do K = 0, kmax
			do M = -K, K
				do N = -K, K
					dr(K,M,N) = reduced(K,M,N) * exp(-ii*(alpha*M+gamma*N))
				enddo
			enddo
		enddo

   end subroutine

!----------------------------------------------------------------
! This function returns the Kronecker delta
!----------------------------------------------------------------
	function delta(A,B)
	integer :: A, B, delta
		delta = 0.d0
		if (A == B) delta = 1.d0
	end function delta

!----------------------------------------------------------------
! This function calculates the t^K_P(i) symbol
!----------------------------------------------------------------
   function tkp(K,P)
	integer :: K, P
	complex(kind=8) :: tkp(0:3), ii

		tkp = 0.d0
		ii = cmplx(0.d0,1.d0)
		if (K == 0) then
			tkp(0) = delta(P,0)
			return
		endif

		if (K == 1) then
			tkp(3) = sqrt(3.d0/2.d0) * delta(P,0)
			return
		endif

		if (K == 2) then
			tkp(0) = sqrt(1.d0/2.d0) * delta(P,0)
			tkp(1) = -sqrt(3.d0)/2.d0 * (delta(P,-2) + delta(P,2))
			tkp(2) = ii * sqrt(3.d0)/2.d0 * (delta(P,-2) - delta(P,2))
			return
		endif
	end function tkp

!----------------------------------------------------------------
! This function calculates the tensor T^K_P(i,Omega) taking into account a double rotation
!----------------------------------------------------------------
   subroutine tkq2(alpha,beta,gamma,alphap,betap,gammap,tr)
	real(kind=8) :: alpha, beta, gamma, alphap, betap, gammap
	complex(kind=8) :: tr(0:2,-2:2,0:3), tkpi(0:3), dr1(0:2,-2:2,-2:2), dr2(0:2,-2:2,-2:2), suma
	integer :: i, K, Q, P, Qp

		call rotmatsu(2,alpha,beta,gamma,dr1)
		call rotmatsu(2,alphap,betap,gammap,dr2)

		tr = 0.d0

		do i = 0, 3
			do K = 0, 2
				do Q = -K, K
					suma = 0.d0
					do P = -K, K
						tkpi = tkp(K,P)
						do Qp = -K, K
							suma = suma + tkpi(i) * dr1(K,P,Qp) * dr2(K,Qp,Q)
						enddo
					enddo
					tr(K,Q,i) = suma
				enddo
			enddo
		enddo

	end subroutine tkq2

!----------------------------------------------------------------
! Calculates part of eq. 7.41 from Landi Degl'Innocenti & Landolfi (2004)
!----------------------------------------------------------------
   subroutine bigpar(l2,is2,k2,kp2,j2,jp2,js2,jt2,x)
	integer :: l2, is2, k2, kp2, j2, jp2, js2, jt2
	real(kind=8) :: x, x1, x2, x3, g

      x = 0.d0
      if (jp2 /= jt2 .and. j2 /= js2) return
      if (jp2 == jt2) then
      	if (iabs(j2-js2) > 2 .or. j2+js2 == 0) then
			else
      		call gamma(l2,is2,j2,js2,g)
      		x1 = g
      		x2 = w6js(k2,kp2,2,js2,j2,jp2)
      		x = x1*x2
			endif
		endif

      if (j2 /= js2) return
      if (iabs(jt2-jp2) > 2 .or. jt2+jp2 == 0) return
      call gamma(l2,is2,jt2,jp2,g)
      x1 = g
      x2 = 1.d0
      if (mod(k2-kp2,4) /= 0) x2 = -1.d0
      x3 = w6js(k2,kp2,2,jt2,jp2,j2)
      x = x + x1*x2*x3
	end subroutine bigpar

!----------------------------------------------------------------
! Calculates eq. 7.42 from Landi Degl'Innocenti & Landolfi (2004)
!----------------------------------------------------------------
   subroutine gamma(l2,is2,j2,jp2,g)
	integer :: l2, is2, j2, jp2
	real(kind=8) :: g
	real(kind=8) :: x1, x2, x3

      g = 0.d0
      if (j2 == jp2) then
			g = g + dsqrt(dfloat(j2*(j2+2)*(j2+1))/4.d0)
		endif
      x1 = 1.d0
      if (mod(2+l2+is2+j2,4) /= 0) x1 = -1.d0
      x2 = dsqrt(dfloat((j2+1)*(jp2+1)*is2*(is2+2)*(is2+1))/4.d0)
      x3 = w6js(j2,jp2,2,is2,is2,l2)
      g = g + x1*x2*x3
	end subroutine gamma

!----------------------------------------------------------------
! Calculates the factorial of the integers up to 301 and save it into the fact array
!----------------------------------------------------------------
   subroutine factrl
	integer :: i
      fact(0) = 1.d0
      do i=1,301
			fact(I) = fact(I-1) * dble(I)
		enddo
   END subroutine factrl

! ---------------------------------------------------------
! LU decomposition of a matrix
!  INPUT:
!		- a is the matrix to decompose
!
!  OUTPUT:
!		- a is the LU decomposition of a
!		- indx is a vector that records the row permutation effected by the partial pivoting
!		- d takes values +1/-1 depending on whether the number of row interchanges was odd or even
! ---------------------------------------------------------
	subroutine ludcmp(a,indx,d)
	integer, INTENT(INOUT) :: indx(:)
	real*8, INTENT(INOUT) :: a(:,:), d
	real*8, parameter :: TINY = 1.d-20
	integer :: i, imax, j, k, n
	real*8 :: aamax, dum, sum, vv(size(a,1))
	integer :: values_start(8), values_end(8)
	real(kind=8) :: start_time, end_time

!		call date_and_time(values=values_start)

		d = 1.d0
		n = size(a,1)

		do i = 1, n
			aamax = 0.d0
			aamax = maxval(dabs(a(i,:)))
			if (aamax == 0.d0) then
				print *, 'Singular matrix in LU decomposition'
				stop
			endif
			vv(i) = 1.d0 / aamax
		enddo

		do j = 1, n
			do i = 1, j-1
				sum = a(i,j)
				do k = 1, i-1
					sum = sum - a(i,k) * a(k,j)
				enddo
				a(i,j) = sum
			enddo
			aamax = 0.d0
			do i = j, n
				sum = a(i,j)
				do k = 1, j-1
					sum = sum - a(i,k) * a(k,j)
				enddo
				a(i,j) = sum
				dum = vv(i) * dabs(sum)
				if (dum >= aamax) then
					imax = i
					aamax = dum
				endif
			enddo
			if (j /= imax) then
				do k = 1, n
					dum = a(imax,k)
					a(imax,k) = a(j,k)
					a(j,k) = dum
				enddo
				d = -d
				vv(imax) = vv(j)
			endif
			indx(j) = imax
			if (a(j,j) == 0.d0) a(j,j) = TINY
			if (j /= n) then
				dum = 1.d0 / a(j,j)
				do i = j+1, n
					a(i,j) = a(i,j) * dum
				enddo
			endif
		enddo

!		call date_and_time(values=values_end)
!		start_time = values_start(5) * 3600 + values_start(6) * 60 + values_start(7) + 0.001 * values_start(8)
!		end_time = values_end(5) * 3600 + values_end(6) * 60 + values_end(7) + 0.001 * values_end(8)
!		print *, 'Time : ', end_time-start_time

	end subroutine ludcmp

! ---------------------------------------------------------
! Solves the set of equations AX=b where A is the LU decomposition of a matrix
!  INPUT:
!		- a is the LU decomposition of the system matrix
!		- b is the right hand side vector of the system
! 		- indx is the vector returned by ludcmp
!  OUTPUT:
!		- b is the solution of the system
! ---------------------------------------------------------
	subroutine lubksb(a,indx,b)
	real*8, INTENT(IN) :: a(:,:)
	real*8, INTENT(INOUT) :: b(:)
	integer, INTENT(IN) :: indx(:)
	integer :: i, ii, n, j, ll
	real*8 :: sum
		n = size(a,1)
		ii = 0
		do i = 1, n
			ll = indx(i)
			sum = b(ll)
			b(ll) = b(i)
			if (ii /= 0) then
				do j = ii, i-1
					sum = sum - a(i,j) * b(j)
				enddo
			else if (sum /= 0.d0) then
				ii = i
			endif
			b(i) = sum
		enddo
		do i = n, 1, -1
			sum = b(i)
			do j = i+1, n
				sum = sum - a(i,j) * b(j)
			enddo
			b(i) = sum / a(i,i)
		enddo
	end subroutine lubksb

! ---------------------------------------------------------
! Read lines that don't mind in a file
! INPUT:
!     - unit : unit to whom the file is associated
!     - n : number of lines to read
! ---------------------------------------------------------
   subroutine lb(unit, n)
   integer, INTENT(IN) :: unit, n
   integer :: i
   		  do i = 1, n
   					 read(unit,*)
   		  enddo
   end subroutine lb

! ---------------------------------------------------------
! Returns the sizes of the eigenvectors of the Paschen Back effect
! ---------------------------------------------------------
	subroutine paschen_size(l2,is2,j2max,njlargest)
	integer :: l2, is2
	integer :: m2, j2max, j2min, j2a, j2b, nj, njlargest

		j2min = abs(l2-is2)
		j2max = l2 + is2

! First find the largest value of nj
		njlargest = 0
		do m2 = -j2max, j2max, 2
			j2a = abs(m2)
			if (j2a < j2min) j2a = j2min
			j2b = j2max
			nj = 1 + (j2b-j2a) / 2
			if (nj > njlargest) njlargest = nj
		enddo
	end subroutine paschen_size

! ---------------------------------------------------------
! Calculates the eigenvalues and eigenvectors of the Paschen-Back matrix
! INPUT:
!     l2 = 2*L
!     is2 = 2*S
!     e0 = a vector containing the energy in cm^-1 of the level J (Jmin <= J <= Jmax) organized in terms of the J2 value
!          ex. e0(3) is the energy of the J=3/2 level. (Some elements remain undefined)
!     b = magnetic field in gauss
! OUTPUT:
!     njlev = number of levels for each m2
!     e = a vector containing the eigenvalues of the matrix. It is organized as e(m2,jsmall), where m2 is 2*M and
!         jsmall is a quantum number which takes values between 1 and the number of J compatible with M
!     c = eigenvectors, organized as c(m2,jsmall,J2)
! ---------------------------------------------------------
	subroutine paschen(l2,is2,e0,B,njlev,e,c,j2m,njl)
	integer :: l2, is2, j2m, njl
	integer :: njlev(-j2m:j2m)
	real(kind=8) :: e(-j2m:j2m,njl), c(-j2m:j2m,njl,0:j2m)
	real(kind=8) :: b, e0(0:jlimit2), bb
	real(kind=8), allocatable :: d(:), ud(:), z(:,:)
	integer :: m2, j2max, j2min, j2a, j2b, nj, j2vero, j, jp, jp2vero, njlargest, i, jsmall, jaux, j2
	real(kind=8) :: x1, x2, x3, x4

		bb = PE / (4.d0*PI*PME*PC**2) * B
		j2min = abs(l2-is2)
		j2max = l2 + is2

		do m2 = -j2max, j2max, 2

! Build the matrix which corresponds to the value of M
			j2a = abs(m2)
			if (j2a < j2min) j2a = j2min
			j2b = j2max
			nj = 1 + (j2b-j2a) / 2
			njlev(m2) = nj

			allocate(d(nj))
			d = 0.d0
			allocate(ud(nj))
			ud = 0.d0
			allocate(z(nj,nj))

! First the diagonal part
			do j = 1, nj
				j2vero = j2a+2*(j-1)
      		d(j) = e0(j2vero)
      		x1 = 1.d0
      		if (mod(2*j2vero+l2+is2+m2,4) /= 0) x1 = -1.d0
      		x2 = dfloat(j2vero+1)*dsqrt(dfloat(is2*(is2+2)*(is2+1))/4.d0)
      		x3 = w3js(j2vero,j2vero,2,-m2,m2,0)
      		x4 = w6js(j2vero,j2vero,2,is2,is2,l2)
      		d(j) = d(j) + bb * (dfloat(m2)/2.d0+x1*x2*x3*x4)
			enddo

! If we use the Pashen-Back regime
			if (use_paschen_back == 1) then
! Then the upper diagonal (only if the matrix is at least 2x2)
				if (nj /= 1) then
					do jp = 2, nj
      				jp2vero = j2a+2*(jp-1)
      				j2vero = jp2vero-2
      				x1 = 1.d0
      				if (mod(j2vero+jp2vero+l2+is2+m2,4) /= 0) x1 = -1.d0
      				x2 = dsqrt(dfloat((j2vero+1)*(jp2vero+1)*is2*(is2+2)*(is2+1))/4.d0)
      				x3 = w3js(j2vero,jp2vero,2,-m2,m2,0)
      				x4 = w6js(j2vero,jp2vero,2,is2,is2,l2)
      				ud(jp) = bb*x1*x2*x3*x4
					enddo
				endif
			endif

! Initializes the matrix z with the eigenvectors
			z = 0.d0
			do i = 1, nj
				z(i,i) = 1.d0
			enddo

! Diagonalize the matrix
			call tqli(d,ud,nj,nj,z)

			do jsmall = 1, nj
				e(m2,jsmall) = d(jsmall)
				do jaux = 1, nj
					j2 = j2a + 2*(jaux-1)
					c(m2,jsmall,j2) = z(jaux,jsmall)
				enddo
			enddo

			deallocate(d)
			deallocate(ud)
			deallocate(z)

		enddo

	end subroutine paschen

! ---------------------------------------------------------
! Diagonalizes a tridiagonal matrix
! ---------------------------------------------------------
	subroutine tqli(d,e,n,np,z)
	integer :: n, np, i, l, m, iter, k
	real(kind=8) :: d(np),e(np),z(np,np), g, r, s, c, p, f, b, dd

		if (n.gt.1) then
      do 11 i=2,n
      e(i-1)=e(i)
 11   continue
      e(n)=0.d0
      do 15 l=1,n
      iter=0
 1    do 12 m=l,n-1
      dd=dabs(d(m))+dabs(d(m+1))
      if(dabs(e(m))+dd.eq.dd) go to 2
 12   continue
      m=n
 2    if(m.ne.l) then
      if(iter.eq.30) then
      	print *, 'too many iterations'
      	stop
      endif
      iter=iter+1
      g=(d(l+1)-d(l))/(2.d0*e(l))
      r=dsqrt(g**2+1.d0)
      g=d(m)-d(l)+e(l)/(g+dsign(r,g))
      s=1.d0
      c=1.d0
      p=0.d0
      do 14 i=m-1,l,-1
      f=s*e(i)
      b=c*e(i)
      if(dabs(f).ge.dabs(g)) then
      c=g/f
      r=dsqrt(c**2+1.d0)
      e(i+1)=f*r
      s=1.d0/r
      c=c*s
      else
      s=f/g
      r=dsqrt(s**2+1.d0)
      e(i+1)=g*r
      c=1.d0/r
      s=s*c
      endif
      g=d(i+1)-p
      r=(d(i)-g)*s+2.d0*c*b
      p=s*r
      d(i+1)=g+p
      g=c*r-b
      do 13 k=1,n
      f=z(k,i+1)
      z(k,i+1)=s*z(k,i)+c*f
      z(k,i)=c*z(k,i)-s*f
 13   continue
 14   continue
      d(l)=d(l)-p
      e(l)=g
      e(m)=0.d0
      go to 1
      endif
 15   continue
      endif
      return
      end subroutine tqli

! ---------------------------------------------------------
! Sets the imaginary part of a complex to a given number
! ---------------------------------------------------------
		subroutine set_imaginary(A,b)
		complex(kind=8) :: A
		real(kind=8) :: b
			A = cmplx(real(A),b)
		end subroutine set_imaginary

! ---------------------------------------------------------
! Sets the imaginary part of a complex to a given number
! ---------------------------------------------------------
		subroutine set_real(A,b)
		complex(kind=8) :: A
		real(kind=8) :: b
			A = cmplx(b,aimag(A))
		end subroutine set_real

!--------------------------------------------------------------
! Generates Voigt and anomalous dispersion profiles using a different damping for each frequency
! point
! See Humlicek (1982) JQSRT 27, 437
!--------------------------------------------------------------
	function profile(da, dv)
	real(kind=8) :: da
	real(kind=8) :: dv(:)
	complex(kind=8) :: profile(size(dv))
	complex(kind=8) :: w4, z, t, u, v4
	real(kind=8) :: s
	integer :: i, n

		n = size(dv)
		if (da < 1.d-4) da = 0.d0
		do i = 1, n
! There is a problem for very small dampings. We force it to be 0 using a threshold
			z = cmplx(dv(i), da)
			t = cmplx(da, -dv(i))
			s = dabs(dv(i)) + da
			u = t*t

			if (s >= 15.d0) then
				w4 = t * 0.5641896d0 / (0.5d0+u)
			elseif (s >= 5.5) then
				w4 = t*(1.410474d0+u*0.5641896d0)/(0.75d0+u*(3.d0+u))
			elseif (da >= 0.195d0*dabs(dv(i))-0.176d0) then
				w4 = (16.4955d0+t*(20.20933d0+t*(11.96482d0+t*(3.778987d0+t*0.5642236d0)))) / &
					(16.4955d0+t*(38.82363d0+t*(39.27121d0+t*(21.69274d0+t*(6.699398d0+t)))))
			else
				w4 = t*(36183.31d0-u*(3321.9905d0-u*(1540.787d0-u*(219.0313d0-u*(35.76683d0-&
					u*(1.320522d0-u*0.56419d0))))))
				v4 = (32066.6d0-u*(24322.84d0-u*(9022.228d0-u*(2186.181d0-u*(364.2191d0-&
					u*(61.57037d0-u*(1.841439d0-u)))))))
				w4 = exp(u) - w4/v4
			endif
			profile(i) = w4
		enddo

	end function profile

! ---------------------------------------------------------
! Given xa(:) and ya(:) which tabulate a function, returns the interpolation using
! linear interpolation of vector x(:) in y(:)
! ---------------------------------------------------------
	subroutine lin_interpol(xa,ya,x,y)
	real*8, INTENT(IN) :: xa(:), ya(:), x(:)
	real*8, INTENT(INOUT) :: y(:)
	integer :: i, n
	integer :: locat(1), loc

		n = size(x)

		do i = 1, n
			locat = minloc(dabs(xa-x(i)))
			loc = locat(1)
			if (loc == 1) then
				y(i) = ya(loc)
			else
				y(i) = (ya(loc)-ya(loc-1))/(xa(loc)-xa(loc-1)) * (x(i)-xa(loc-1)) + ya(loc-1)
			endif
		enddo

	end subroutine lin_interpol

! ---------------------------------------------------------
! Given etaI, etaQ, etaU, etaV, rhoQ, rhoU and rhoV, fill the absorption matrix
! ---------------------------------------------------------
	subroutine fill_absorption_matrix(matrix,etaI,etaQ,etaU,etaV,rhoQ,rhoU,rhoV)
	real(kind=8) :: matrix(4,4), etaI, etaQ, etaU, etaV, rhoQ, rhoU, rhoV

		matrix(1,1) = etaI
		matrix(2,2) = etaI
		matrix(3,3) = etaI
		matrix(4,4) = etaI
		matrix(1,2) = etaQ
		matrix(2,1) = etaQ
		matrix(1,3) = etaU
		matrix(3,1) = etaU
		matrix(1,4) = etaV
		matrix(4,1) = etaV
		matrix(2,3) = rhoV
		matrix(3,2) = -rhoV
		matrix(2,4) = -rhoU
		matrix(4,2) = rhoU
		matrix(3,4) = rhoQ
		matrix(4,3) = -rhoQ

	end subroutine fill_absorption_matrix

!--------------------------------------------------------------
! Inversion of a 4x4 matrix
!--------------------------------------------------------------
	subroutine invert(a)
	real*8 :: a(4,4)
	real*8 :: b(4,4), det, maxim, fabsmax
	! First some tests of singularity
		b = dabs(a)
		maxim = maxval(b)
		fabsmax = 1.d0 / maxim
		if (maxim == 0.d0) then
			print *, 'Singularity in the inversion'
!			stop
		endif

		a = a * fabsmax

   	b(1,1) = a(2,2) * a(3,3) * a(4,4) + a(2,3) * a(3,4) * a(4,2)&
      	+ a(2,4) * a(3,2) * a(4,3) - a(2,2) * a(3,4) * a(4,3)&
	   	- a(2,3) * a(3,2) * a(4,4) - a(2,4) * a(3,3) * a(4,2)
   	b(2,1) = a(2,3) * a(3,1) * a(4,4) + a(2,4) * a(3,3) * a(4,1)&
	   	+ a(2,1) * a(3,4) * a(4,3) - a(2,3) * a(3,4) * a(4,1)&
   		- a(2,4) * a(3,1) * a(4,3) - a(2,1) * a(3,3) * a(4,4)
   	b(3,1) = a(2,4) * a(3,1) * a(4,2) + a(2,1) * a(3,2) * a(4,4)&
	   	+ a(2,2) * a(3,4) * a(4,1) - a(2,4) * a(3,2) * a(4,1)&
   		- a(2,1) * a(3,4) * a(4,2) - a(2,2) * a(3,1) * a(4,4)
   	b(4,1) = a(2,1) * a(3,3) * a(4,2) + a(2,2) * a(3,1) * a(4,3)&
	   	+ a(2,3) * a(3,2) * a(4,1) - a(2,1) * a(3,2) * a(4,3)&
   		- a(2,2) * a(3,3) * a(4,1) - a(2,3) * a(3,1) * a(4,2)
   	b(1,2) = a(3,2) * a(4,4) * a(1,3) + a(3,3) * a(4,2) * a(1,4)&
	   	+ a(3,4) * a(4,3) * a(1,2) - a(3,2) * a(4,3) * a(1,4)&
   		- a(3,3) * a(4,4) * a(1,2) - a(3,4) * a(4,2) * a(1,3)
   	b(2,2) = a(3,3) * a(4,4) * a(1,1) + a(3,4) * a(4,1) * a(1,3)&
	   	+ a(3,1) * a(4,3) * a(1,4) - a(3,3) * a(4,1) * a(1,4)&
   		- a(3,4) * a(4,3) * a(1,1) - a(3,1) * a(4,4) * a(1,3)
   	b(3,2) = a(3,4) * a(4,2) * a(1,1) + a(3,1) * a(4,4) * a(1,2)&
	   	+ a(3,2) * a(4,1) * a(1,4) - a(3,4) * a(4,1) * a(1,2)&
   		- a(3,1) * a(4,2) * a(1,4) - a(3,2) * a(4,4) * a(1,1)
   	b(4,2) = a(3,1) * a(4,2) * a(1,3) + a(3,2) * a(4,3) * a(1,1)&
	   	+ a(3,3) * a(4,1) * a(1,2) - a(3,1) * a(4,3) * a(1,2)&
   		- a(3,2) * a(4,1) * a(1,3) - a(3,3) * a(4,2) * a(1,1)
   	b(1,3) = a(4,2) * a(1,3) * a(2,4) + a(4,3) * a(1,4) * a(2,2)&
   		+ a(4,4) * a(1,2) * a(2,3) - a(4,2) * a(1,4) * a(2,3)&
	   	- a(4,3) * a(1,2) * a(2,4) - a(4,4) * a(1,3) * a(2,2)
   	b(2,3) = a(4,3) * a(1,1) * a(2,4) + a(4,4) * a(1,3) * a(2,1)&
	   	+ a(4,1) * a(1,4) * a(2,3) - a(4,3) * a(1,4) * a(2,1)&
   		- a(4,4) * a(1,1) * a(2,3) - a(4,1) * a(1,3) * a(2,4)
   	b(3,3) = a(4,4) * a(1,1) * a(2,2) + a(4,1) * a(1,2) * a(2,4)&
	   	+ a(4,2) * a(1,4) * a(2,1) - a(4,4) * a(1,2) * a(2,1)&
   		- a(4,1) * a(1,4) * a(2,2) - a(4,2) * a(1,1) * a(2,4)
   	b(4,3) = a(4,1) * a(1,3) * a(2,2) + a(4,2) * a(1,1) * a(2,3)&
	   	+ a(4,3) * a(1,2) * a(2,1) - a(4,1) * a(1,2) * a(2,3)&
   		- a(4,2) * a(1,3) * a(2,1) - a(4,3) * a(1,1) * a(2,2)
   	b(1,4) = a(1,2) * a(2,4) * a(3,3) + a(1,3) * a(2,2) * a(3,4)&
	   	+ a(1,4) * a(2,3) * a(3,2) - a(1,2) * a(2,3) * a(3,4)&
   		- a(1,3) * a(2,4) * a(3,2) - a(1,4) * a(2,2) * a(3,3)
   	b(2,4) = a(1,3) * a(2,4) * a(3,1) + a(1,4) * a(2,1) * a(3,3)&
	   	+ a(1,1) * a(2,3) * a(3,4) - a(1,3) * a(2,1) * a(3,4)&
   		- a(1,4) * a(2,3) * a(3,1) - a(1,1) * a(2,4) * a(3,3)
   	b(3,4) = a(1,4) * a(2,2) * a(3,1) + a(1,1) * a(2,4) * a(3,2)&
	   	+ a(1,2) * a(2,1) * a(3,4) - a(1,4) * a(2,1) * a(3,2)&
   		- a(1,1) * a(2,2) * a(3,4) - a(1,2) * a(2,4) * a(3,1)
   	b(4,4) = a(1,1) * a(2,2) * a(3,3) + a(1,2) * a(2,3) * a(3,1)&
	   	+ a(1,3) * a(2,1) * a(3,2) - a(1,1) * a(2,3) * a(3,2)&
   		- a(1,2) * a(2,1) * a(3,3) - a(1,3) * a(2,2) * a(3,1)

		det = a(1,1) * b(1,1) + a(1,2) * b(2,1) + a(1,3) * b(3,1) + a(1,4) * b(4,1)

		a = b * (fabsmax / det)

	end subroutine invert

!--------------------------------------------------------------
! Returns in b the evolution operator if the propagation matrix is constant
! given by exp(-A*t). It uses Eqs. A5.20 of Landi Degl'Innocenti & Landolfi 2004
!--------------------------------------------------------------
	subroutine evol_operator(a,t,b)
	real(kind=8) :: a(4,4), t
	real(kind=8) :: b(4,4), det, maxim, fabsmax, w(3), y(3), w_y_scalar, sigma
	real(kind=8) :: Lambda_plus, Lambda_minus, tt, Theta, w_mod_sq, y_mod_sq, factor1, factor2, factor3, factor4, w0, temp1, temp2
	real(kind=8), dimension(4,4) :: M1, M2, M3, M4
	integer :: i

		w0 = a(1,1)

		w(1) = a(1,2)
		w(2) = a(1,3)
		w(3) = a(1,4)
		w_mod_sq = sum(w*w)

		y(1) = a(3,4)
		y(2) = a(4,2)
		y(3) = a(2,3)
		y_mod_sq = sum(y*y)

		w_y_scalar = sum(w*y)
		sigma = 1.d0
		if (w_y_scalar < 0.d0) sigma = -1.d0

! Eigenvalues
		temp1 = sqrt((w_mod_sq - y_mod_sq)**2 / 4.d0 + w_y_scalar**2)
		temp2 = (w_mod_sq - y_mod_sq) / 2.d0
		if (temp1+temp2 < 0.d0) then
			Lambda_plus = 0.d0
		else
			Lambda_plus = sqrt(temp1 + temp2)
		endif
		if (temp1-temp2 < 0.d0) then
			Lambda_minus = 0.d0
		else
			Lambda_minus = sqrt(temp1 - temp2)
		endif

		Theta = 2.d0 * temp1
		tt = (w_mod_sq + y_mod_sq) / 2.d0

		M1 = 0.d0
		M2 = 0.d0
		M3 = 0.d0
		M4 = 0.d0

! M1
		do i = 1, 4
			M1(i,i) = 1.d0
		enddo

! M2
		M2(1,2) = Lambda_minus * w(1) - sigma * Lambda_plus * y(1)
		M2(1,3) = Lambda_minus * w(2) - sigma * Lambda_plus * y(2)
		M2(1,4) = Lambda_minus * w(3) - sigma * Lambda_plus * y(3)
		M2(2,1) = M2(1,2)
		M2(2,3) = sigma * Lambda_plus * w(3) + Lambda_minus * y(3)
		M2(2,4) = -sigma * Lambda_plus * w(2) - Lambda_minus * y(2)
		M2(3,1) = M2(1,3)
		M2(3,2) = -M2(2,3)
		M2(3,4) = sigma * Lambda_plus * w(1) + Lambda_minus * y(1)
		M2(4,1) = M2(1,4)
		M2(4,2) = -M2(2,4)
		M2(4,3) = -M2(3,4)

! M3
		M3(1,2) = Lambda_plus * w(1) + sigma * Lambda_minus * y(1)
		M3(1,3) = Lambda_plus * w(2) + sigma * Lambda_minus * y(2)
		M3(1,4) = Lambda_plus * w(3) + sigma * Lambda_minus * y(3)
		M3(2,1) = M3(1,2)
		M3(2,3) = -sigma * Lambda_minus * w(3) + Lambda_plus * y(3)
		M3(2,4) = sigma * Lambda_minus * w(2) - Lambda_plus * y(2)
		M3(3,1) = M3(1,3)
		M3(3,2) = -M3(2,3)
		M3(3,4) = -sigma * Lambda_minus * w(1) + Lambda_plus * y(1)
		M3(4,1) = M3(1,4)
		M3(4,2) = -M3(2,4)
		M3(4,3) = -M3(3,4)

! M4
		M4(1,1) = tt
		M4(1,2) = w(3)*y(2) - w(2)*y(3)
		M4(1,3) = w(1)*y(3) - w(3)*y(1)
		M4(1,4) = w(2)*y(1) - w(1)*y(2)
		M4(2,1) = -M4(1,2)
		M4(2,2) = w(1)**2 + y(1)**2 - tt
		M4(2,3) = w(1)*w(2) + y(1)*y(2)
		M4(2,4) = w(3)*w(1) + y(3)*y(1)
		M4(3,1) = -M4(1,3)
		M4(3,2) = M4(2,3)
		M4(3,3) = w(2)**2 + y(2)**2 - tt
		M4(3,4) = w(2)*w(3) + y(2)*y(3)
		M4(4,1) = -M4(1,4)
		M4(4,2) = M4(2,4)
		M4(4,3) = M4(3,4)
		M4(4,4) = w(3)**2 + y(3)**2 - tt

		M2 = M2 / Theta
		M3 = M3 / Theta
		M4 = 2.d0 * M4 / Theta

		factor1 = 0.5d0*(cosh(t*Lambda_plus) + cos(t*Lambda_minus))
		factor2 = sin(t*Lambda_minus)
		factor3 = sinh(t*Lambda_plus)
		factor4 = 0.5d0*(cosh(t*Lambda_plus) - cos(t*Lambda_minus))

		b = factor1*M1 - factor2*M2 - factor3*M3 + factor4*M4
		b = exp(-t*w0) * b

	end subroutine evol_operator

! ---------------------------------------------------------
! Return the input array but shifted by a sub-pixel amount
! Note that it works only for any length of series which has no prime factor greater than 23
! ---------------------------------------------------------
		function fft_shift(x, sh)
		real(kind=8) :: x(:), fft_shift(size(x)), sh
		complex(kind=8) :: k(size(x)), ff(size(x))
		integer :: nx, i, n21

			nx = size(x)

			do i = 1, nx
         	k(i) = i-1
         enddo
         n21 = nx/2+1
         do i = n21+1, nx
         	k(i) = n21-nx+(i-1-n21)
         enddo

         k = 2.d0*PI*k/(1.d0*nx)
         k = k * cmplx(0.d0,1.d0)

         ff = x

         ff = fft(ff)

         fft_shift = real(fft(ff * exp(k*sh), inv=.TRUE.))

		end function fft_shift

end module maths

max_stars_repo_name: aasensio/pyiacsun

max_issues_repo_name: aasensio/pyiacsun

max_forks_repo_name: aasensio/pyiacsun


text:
! Copyright (C) 2021 Atsushi Togo
! All rights reserved.

! This file is part of kspclib.

! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:

! * Redistributions of source code must retain the above copyright
!   notice, this list of conditions and the following disclaimer.

! * Redistributions in binary form must reproduce the above copyright
!   notice, this list of conditions and the following disclaimer in
!   the documentation and/or other materials provided with the
!   distribution.

! * Neither the name of the kspclib project nor the names of its
!   contributors may be used to endorse or promote products derived
!   from this software without specific prior written permission.

! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.

module kspclib_f08

  use iso_c_binding, only: c_char, c_long, c_double

  implicit none

  private

  interface

     function ksp_get_major_version() bind(c)
       import c_long
       integer(c_long) :: ksp_get_major_version
     end function ksp_get_major_version

     function ksp_get_minor_version() bind(c)
       import c_long
       integer(c_long) :: ksp_get_minor_version
     end function ksp_get_minor_version

     function ksp_get_micro_version() bind(c)
       import c_long
       integer(c_long) :: ksp_get_micro_version
     end function ksp_get_micro_version

     subroutine ksp_get_all_grid_addresses(grid_address, mesh) bind(c)
       import c_long
       integer(c_long), intent(in) :: mesh(3)
       integer(c_long), intent(inout) :: grid_address(3, *)
     end subroutine ksp_get_all_grid_addresses

     subroutine ksp_get_double_grid_address(address_double, address, mesh, &
          is_shift) bind(c)
       import c_long
       integer(c_long), intent(inout) :: address_double(3)
       integer(c_long), intent(in) :: address(3)
       integer(c_long), intent(in) :: mesh(3)
       integer(c_long), intent(in) :: is_shift(3)
     end subroutine ksp_get_double_grid_address

     function ksp_get_double_grid_index(address_double, mesh) bind(c)
       import c_long
       integer(c_long), intent(in) :: mesh(3)
       integer(c_long), intent(in) :: address_double(3)
       integer(c_long) :: ksp_get_double_grid_index
     end function ksp_get_double_grid_index

     subroutine ksp_get_thm_relative_grid_addresses(relative_grid_addresses, &
          rec_lattice) bind(c)
       import c_long, c_double
       integer(c_long), intent(inout) :: relative_grid_addresses(3, 4, 24)
       real(c_double), intent(in) :: rec_lattice(3, 3)
     end subroutine ksp_get_thm_relative_grid_addresses

     function ksp_get_thm_integration_weight(omega, &
          tetrahedra_omegas, function_char) bind(c)
       import c_char, c_double
       real(c_double), intent(in), value :: omega
       real(c_double), intent(in) :: tetrahedra_omegas(4, 24)
       character(kind=c_char), intent(in), value :: function_char
       real(c_double) :: ksp_get_thm_integration_weight
     end function ksp_get_thm_integration_weight

     function ksp_snf_transform_rotations(transformed_rots, &
          rotations, num_rot, D_diag, Q) bind(c)
       import c_long
       integer(c_long), intent(inout) :: transformed_rots(3, 3, *)
       integer(c_long), intent(in) :: rotations(3, 3, *)
       integer(c_long), intent(in), value :: num_rot
       integer(c_long), intent(in) :: D_diag(3)
       integer(c_long), intent(in) :: Q(3, 3)
     end function ksp_snf_transform_rotations

     function ksp_get_snf3x3(D_diag, P, Q, A) bind(c)
       import c_long
       integer(c_long), intent(inout) :: D_diag(3)
       integer(c_long), intent(inout) :: P(3, 3)
       integer(c_long), intent(inout) :: Q(3, 3)
       integer(c_long), intent(in) :: A(3, 3)
       integer(c_long) :: ksp_get_snf3x3
     end function ksp_get_snf3x3

     subroutine ksp_get_all_grgrid_addresses(grid_address, D_diag) bind(c)
       import c_long
       integer(c_long), intent(inout) :: grid_address(3, *)
       integer(c_long), intent(in) :: D_diag(3)
     end subroutine ksp_get_all_grgrid_addresses

     subroutine ksp_get_double_grgrid_address(address_double, address, &
          D_diag, PS) bind(c)
       import c_long
       integer(c_long), intent(inout) :: address_double(3)
       integer(c_long), intent(in) :: D_diag(3)
       integer(c_long), intent(in) :: address(3)
       integer(c_long), intent(in) :: PS(3)
     end subroutine ksp_get_double_grgrid_address

     function ksp_get_grgrid_index(address, D_diag) bind(c)
       import c_long
       integer(c_long), intent(in) :: address(3)
       integer(c_long), intent(in) :: D_diag(3)
       integer(c_long) :: ksp_get_grgrid_index
     end function ksp_get_grgrid_index

     function ksp_get_double_grgrid_index(address_double, D_diag, PS) bind(c)
       import c_long
       integer(c_long), intent(in) :: address_double(3)
       integer(c_long), intent(in) :: D_diag(3)
       integer(c_long), intent(in) :: PS(3)
       integer(c_long) :: ksp_get_double_grgrid_index
     end function ksp_get_double_grgrid_index

     subroutine ksp_get_grgrid_address_from_index(address, grid_index, D_diag) &
          bind(c)
       import c_long
       integer(c_long), intent(inout) :: address(3)
       integer(c_long), intent(in), value :: grid_index
       integer(c_long), intent(in) :: D_diag(3)
     end subroutine ksp_get_grgrid_address_from_index

     function ksp_rotate_grgrid_index(grid_index, rotation, D_diag, PS) bind(c)
       import c_long
       integer(c_long), intent(in), value :: grid_index
       integer(c_long), intent(in) :: rotation(3, 3)
       integer(c_long), intent(in) :: D_diag(3)
       integer(c_long), intent(in) :: PS(3)
       integer(c_long) :: ksp_rotate_grgrid_index
     end function ksp_rotate_grgrid_index

     subroutine ksp_get_ir_grgrid_map(ir_grid_indices, rotations, num_rot, &
          D_diag, PS) bind(c)
       import c_long
       integer(c_long), intent(inout) :: ir_grid_indices(*)
       integer(c_long), intent(in) :: rotations(3, 3, *)
       integer(c_long), intent(in), value :: num_rot
       integer(c_long), intent(in) :: D_diag(3)
       integer(c_long), intent(in) :: PS(3)
     end subroutine ksp_get_ir_grgrid_map

     function ksp_get_reciprocal_point_group(rec_rotations, rotations, &
          num_rot, is_time_reversal) bind(c)
       import c_long
       integer(c_long), intent(inout) :: rec_rotations(3, 3, 48)
       integer(c_long), intent(in) :: rotations(3, 3, *)
       integer(c_long), intent(in), value :: num_rot
       integer(c_long), intent(in), value :: is_time_reversal
       integer(c_long) :: ksp_get_reciprocal_point_group
     end function ksp_get_reciprocal_point_group
  end interface

  public :: ksp_get_major_version,  ksp_get_minor_version, &
       ksp_get_micro_version, ksp_get_all_grid_addresses, ksp_get_snf3x3, &
       ksp_get_double_grid_address, ksp_get_double_grid_index, &
       ksp_get_thm_relative_grid_addresses, ksp_get_thm_integration_weight, &
       ksp_snf_transform_rotations, ksp_get_all_grgrid_addresses, &
       ksp_get_double_grgrid_address, ksp_get_grgrid_index, &
       ksp_get_double_grgrid_index, ksp_get_grgrid_address_from_index, &
       ksp_rotate_grgrid_index, ksp_get_ir_grgrid_map, &
       ksp_get_reciprocal_point_group

end module kspclib_f08

max_stars_repo_name: atztogo/kspclib

max_issues_repo_name: atztogo/kspclib

max_forks_repo_name: atztogo/kspclib


text:
module TMPL_Demo_M

use ISO_C_BINDING

implicit none

private

interface
    function TMPL_c_function(j) result(i) BIND(C,name="TMPL_c_function")
    use ISO_C_BINDING
    integer(C_INT), value :: j
    integer(C_INT) :: i
    end function
end interface

public :: TMPL_bound_function

contains

function TMPL_bound_function(j) result(i)
    integer, intent(in) :: j
    integer :: i
    write(*,*)"here with j=",j
    i=TMPL_c_function(int(j,C_INT))
    write(*,*)'after'
end function

end module

max_stars_repo_name: wawiesel/BOTG-SkeletonRepo

max_issues_repo_name: wawiesel/BOTG-SkeletonRepo

max_forks_repo_name: wawiesel/BOTG-SkeletonRepo


program neqDMFT
  USE NEQ_DMFT_IPT
  USE SCIFOR
  USE DMFT_TOOLS
  implicit none
  integer                               :: i,j,ik,itime,iloop,ix,iy,Lk
  logical                               :: converged
  real(8)                               :: wband
  character(len=16)                     :: finput
  type(kb_contour_gf)                   :: Sbath
  type(kb_contour_gf)                   :: Gloc
  type(kb_contour_gf)                   :: Gwf
  type(kb_contour_sigma)                :: Sigma
  !
  type(kb_contour_gf)                   :: Ker
  type(kb_contour_gf),allocatable       :: Gk(:)
  type(kb_contour_dgf),allocatable      :: dGk(:),dGk_old(:)
  complex(8),dimension(:,:),allocatable :: Ham
  complex(8)                            :: sigma_gtr
  !RESULTS:
  real(8),dimension(:),allocatable      :: epsik,wt
  real(8),dimension(:,:),allocatable    :: nk
  real(8)                               :: D,de,intwt

  !READ THE INPUT FILE (in vars_global):
  call parse_cmd_variable(finput,"FINPUT",default='inputNEQ.conf')
  call parse_input_variable(wband,"wband",finput,default=2d0,comment="half-bandwidth W=2t")
  call parse_input_variable(Lk,"Lk",finput,default=100,comment="Number of energy levels for Bethe DOS")
  call read_input_init(trim(finput))

  !BUILD TIME GRIDS AND NEQ-PARAMETERS:
  call allocate_kb_contour_params(cc_params,Ntime,Ntau,Niw)
  call setup_kb_contour_params(cc_params,dt,beta)

  !BUILD THE LATTICE STRUCTURE (in lib/square_lattice):
  allocate(epsik(Lk),wt(Lk))
  epsik = linspace(-wband,wband,Lk,mesh=de)
  do i=1,Lk
     wt(i) = dens_bethe(epsik(i),wband)
  enddo
  call splot("DOSbethe.ipt",epsik,wt)
  wt=wt/sum(wt)

  !ALLOCATE ALL THE FUNCTIONS INVOLVED IN THE CALCULATION:
  call allocate_kb_contour_sigma(Sigma,cc_params) !Self-Energy function
  call allocate_kb_contour_gf(Gloc,cc_params)     !Local Green's function
  call allocate_kb_contour_gf(Gwf,cc_params)      !Local Weiss-Field function
  call allocate_kb_contour_gf(Ker,cc_params)
  allocate(Gk(Lk),dGk(Lk),dGk_old(Lk))
  do ik=1,Lk
     call allocate_kb_contour_gf(Gk(ik),cc_params)
     call allocate_kb_contour_dgf(dGk(ik),cc_params)
     call allocate_kb_contour_dgf(dGk_old(ik),cc_params)
  end do
  allocate(ham(cc_params%Ntime,Lk))
  allocate(nk(cc_params%Ntime,Lk))

  !READ OR GUESS THE INITIAL WEISS FIELD G0 (t=t'=0)
  cc_params%Nt=1
  Gloc = zero
  call neq_continue_equilibirum(Gwf,Gk,dGk,Gloc,Sigma,epsik,wt,cc_params)
  do ik=1,Lk
     nk(1,ik)=dimag(Gk(ik)%less(1,1))
  enddo
  do i=1,cc_params%Ntime
     do ik=1,Lk
        ham(i,ik)=epsik(ik)
     enddo
  enddo
  call measure_observables(Gloc,Sigma,cc_params)

  !START THE TIME_STEP LOOP  1<t<=Nt
  !AT EACH TIME_STEP PERFORM A FULL DMFT CALCULATION:
  do itime=2,cc_params%Ntime
     print*,""
     print*,"time step=",itime
     cc_params%Nt=itime

     !prepare the weiss-field at this actual time_step for DMFT:
     call  extrapolate_kb_contour_gf(Gwf,cc_params)
     call  extrapolate_kb_contour_gf(Gloc,cc_params)
     do ik=1,Lk
        dGk_old(ik) = dGk(ik)
     enddo

     !DMFT LOOP:
     iloop=0;converged=.false.
     do while(.not.converged.AND.iloop<nloop)
        iloop=iloop+1
        write(*,"(A,I4,A1)",advance='no')"dmft loop=",iloop," "

        !IMPURITY SOLVER: IPT.
        call neq_solve_ipt(Gwf,Gloc,Sigma,cc_params)

        !PERFORM THE SELF_CONSISTENCY: get local GF + update Weiss-Field
        do ik=1,Lk
           call vide_kb_contour_gf(Ham(:,ik),Sigma,Gk(ik),dGk_old(ik),dGk(ik),cc_params)
        enddo
        call sum_kb_contour_gf(Gk(:),wt(:),Gloc,cc_params)

        !update the weiss field by solving the integral equation:
        ! G0  = Q + K*G0 , with K = -G*\Sigma and Q = G
        call convolute_kb_contour_gf(Gloc,Sigma,Ker,cc_params,dcoeff=-1.d0)
        call vie_kb_contour_gf(Gloc,Ker,Gwf,cc_params)
        !
        !CHECK CONVERGENCE
        converged = convergence_check(Gwf,cc_params)
     enddo

     !EVALUATE AND PRINT THE RESULTS OF THE CALCULATION
     call measure_observables(Gloc,Sigma,cc_params)
     forall(ik=1:Lk)nk(itime,ik)=dimag(Gk(ik)%less(itime,itime))
  enddo

  !EVALUATE AND PRINT OTHER RESULTS OF THE CALCULATION
  call splot3d("nkVSepsikVStime.neqipt",cc_params%t,epsik,nk)
  call plot_kb_contour_sigma("Sigma",Sigma,cc_params)
  call plot_kb_contour_gf("Gloc",Gloc,cc_params)
  call plot_kb_contour_gf("G0",Gwf,cc_params)
  print*,"BRAVO"

contains

  !+-------------------------------------------------------------------+
  !PURPOSE  : check convergence of the calculation:
  !+-------------------------------------------------------------------+
  function convergence_check(G,params) result(converged)
    type(kb_contour_gf)                 :: G
    type(kb_contour_params)             :: params
    logical                             :: converged
    complex(8),dimension(:),allocatable :: test_func
    integer                             :: i,N,L,Ntot
    !
    N   = params%Nt                 !<== work with the ACTUAL size of the contour
    L   = params%Ntau
    !
    Ntot=2*N+L+1
    allocate(test_func(Ntot))
    test_func=zero
    do i=1,N
       test_func(i)  = G%ret(N,i)
       test_func(N+i)= G%less(N,i)
    enddo
    do i=0,L
       test_func(2*N+i+1)=G%lmix(N,i)
    enddo
    converged=check_convergence(test_func,dmft_error,Nsuccess,nloop)
    deallocate(test_func)
    !if(isnan(err))stop "Aborted convergence: error=NaN"
  end function convergence_check

end PROGRAM neqDMFT


text:
program example_mnist

  ! A training example with the MNIST dataset.
  ! Uses stochastic gradient descent and mini-batch size of 100.
  ! Can be run in serial or parallel mode without modifications.

  use mod_kinds, only: ik, rk
  use mod_mnist, only: label_digits, load_mnist
  use mod_network, only: network_type

  implicit none

  real(rk), allocatable :: tr_images(:,:), tr_labels(:)
  real(rk), allocatable :: te_images(:,:), te_labels(:)

  type(network_type) :: net

  integer(ik) :: i, n, num_epochs
  integer(ik) :: batch_size

  call load_mnist(tr_images, tr_labels, te_images, te_labels)

  net = network_type([size(tr_images, dim=1), 10, size(label_digits(tr_labels), dim=1)])

  batch_size = 100
  num_epochs = 10

  if (this_image() == 1) print '(a,f5.2,a)', 'Initial accuracy: ', &
    net % accuracy(te_images, label_digits(te_labels)) * 100, ' %'

  call net % train(tr_images, label_digits(tr_labels), 3._rk, num_epochs, batch_size)

  if (this_image() == 1) print '(a,f5.2,a)', 'Epochs done, Accuracy: ', &
    net % accuracy(te_images, label_digits(te_labels)) * 100, ' %'

end program example_mnist

max_stars_repo_name: GridIO/neural-fortran

max_issues_repo_name: GridIO/neural-fortran

max_forks_repo_name: GridIO/neural-fortran


!=============================================================================!
        program conv
!
!  Convert Fortran unformatted double PLOT3D file to a binary float
!  PLOT3D file on the SGI
!
!  Default conversion on the qfile (or use -q) -- To translate a grid file
!  use -g option
!
!  Author:   S. Scott Collis
!
!  Date:     1-4-97
!
!=============================================================================!
        use const
        implicit none

!.... flow data

        real, allocatable :: q(:,:,:,:)

!.... mesh

        real, allocatable :: x(:,:,:), y(:,:,:), z(:,:,:)

!.... parameters

        integer :: lstep, nx, ny, nz, ndof=5
        real  :: Ma, Re, Pr, time

!.... local variables

        integer :: i, j, k, idof

        character*80 file, temp
        integer :: iloc, iend
#ifndef __GFORTRAN__
        integer, external :: iargc
#endif
        real :: tmp

        logical :: grid=.false., qfile=.true., twod=.false.

        integer, parameter :: mfile = 2
        integer :: narg, iarg, nfile=0, ifile(mfile)
        character*80 :: arg
!=============================================================================!
!.... parse the argument list

        narg = iargc()
        do iarg = 1, narg
          call getarg(iarg,arg)
          if (arg(1:1) .ne. '-') then
            nfile = nfile + 1
            if (nfile .gt. mfile) then
              write(*,*) '>> Error in argument list, too many file names'
              call exit(1)
            end if
            ifile(nfile) = iarg
          else
            select case (arg(1:2))
            case ('-g')
              grid = .true.
              qfile = .false.
            case ('-2')
              twod = .true.
            case ('-h')
              write(*,"('-----------------------------------------------')")
              write(*,"('Usage:  conv [options] [file1] [file2]')")
              write(*,"('-----------------------------------------------')")
              write(*,"(' Convert an unformated double precision PLOT3D')")
              write(*,"(' file to a binary single precision PLOT3D file')")
              write(*,"('-----------------------------------------------')")
              write(*,"('   -h:  this help')")
              write(*,"('   -g:  convert a grid file')")
              write(*,"('   -2:  make a 2-D Plot3d file')")
!             write(*,"('   -q:  convert q data file')")
              write(*,"('-----------------------------------------------')")
              call exit(0)
            case default
              write(*,"('Argument ',i2,' ignored.')") iarg
            end select
          end if
        end do

!.... read in the unformatted PLOT3D grid file

        if (grid) then

        if ( nfile .gt. 0 ) then
          call getarg(ifile(1),temp)
          file = temp
        else
          file = 'grid.dat'
          write (*,"('Enter unformatted grid file name [',a,']? ',$)") &
            file(1:index(file,' ')-1)
          read (*,"(a80)") temp
          if (temp(1:1) .ne. ' ') file = temp
        end if
 10     open(unit=10, file=file, form='unformatted', status='old', err=20)
        goto 30
 20     write (*,"('>> Error opening [',a,'] ',/)") file(1:index(file,' '))
        write (*,"('Enter file name [',a,']? ',$)") file(1:index(file,' '))
        read (*,"(a80)") temp
        if (temp(1:1) .ne. ' ') file = temp
        goto 10
 30     continue

        open(unit=10,file=file,form='unformatted',status='old')
        read(10) nx, ny, nz
!       write(*,*) nx,ny,nz
        allocate( x(nx,ny,nz), y(nx,ny,nz), z(nx,ny,nz) )
        read(10) (((x(i,j,k), i = 1, nx), j = 1, ny), k = 1, nz), &
                 (((y(i,j,k), i = 1, nx), j = 1, ny), k = 1, nz), &
                 (((z(i,j,k), i = 1, nx), j = 1, ny), k = 1, nz)
        close(10)

!.... write out the binary PLOT3D grid file

        iloc = index(file,'.dat')
        iend = index(file,' ')-1
        if (iloc .ne. 0) then
          temp = file(1:iloc)//'bin'
        else
          temp = file(1:iend)//'.bin'
        end if
        if ( nfile .eq. 0 ) then
          file = temp
          write (*,"(/,'Enter binary PLOT3D file name [',a,']? ',$)")   &
                  file(1:index(file,' ')-1)
          read (*,"(a80)") temp
          if (temp(1:1) .ne. ' ') file = temp
        else if ( nfile .eq. 1 ) then
          file = temp
        else if ( nfile .eq. 2 ) then
          call getarg(ifile(2),temp)
          file = temp
        end if
        iloc = index(file,' ')
        file = file(1:(iloc-1))//char(0)

        if (twod) then
          call wgrid2d(x, y, nx, ny, nx, file)
        else
          call wgrid  (x, y, z, nx, ny, nz, nx, ny, file)
        end if

        end if

!.... read the unformatted Plot3d q-file

        if (qfile) then

        if ( nfile .gt. 0 ) then
          call getarg(ifile(1),temp)
          file = temp
        else
          file = 'output.q.1'
          write (*,"('Enter unformatted qfile name [',a,']? ',$)") &
            file(1:index(file,' '))
          read (*,"(a80)") temp
          if (temp(1:1) .ne. ' ') file = temp
        end if
 50     open(unit=10, file=file, form='unformatted', status='old', err=60)
        goto 70
 60     write (*,"('>> Error opening [',a,'] ',/)") file(1:index(file,' '))
        write (*,"('Enter file name [',a,']? ',$)") file(1:index(file,' '))
        read (*,"(a80)") temp
        if (temp(1:1) .ne. ' ') file = temp
        goto 50
 70     continue

        read(10) nx, ny, nz
        read(10) Ma, Pr, Re, time
        allocate ( q(nx,ny,nz,ndof) )
        read(10) ((((q(i,j,k,idof), i = 1, nx), j = 1, ny),  &
                                    k = 1, nz), idof = 1, ndof)
        close(10)

!.... write out the binary Plot3d q-file

        iloc = index(file,'.dat')
        iend = index(file,' ')-1
        if (iloc .ne. 0) then
          temp = file(1:iloc)//'bin'
        else
          temp = file(1:iend)//'.bin'
        end if

        iloc = index(file,'.q.')
        iend = index(file,' ')-1
        if (iloc .ne. 0) then
          temp = file(1:iloc)//'b'//file(iloc+2:iend)
        else
          iloc = index(file,'.dat')
          iend = index(file,' ')-1
          if (iloc .ne. 0) then
            temp = file(1:iloc)//'bin'
          else
            temp = file(1:iend)//'.bin'
          end if
        end if
        if ( nfile .eq. 0 ) then
          file = temp
          write (*,"(/,'Enter binary PLOT3D qfile name [',a,']? ',$)")  &
                  file(1:index(file,' '))
          read (*,"(a80)") temp
          if (temp(1:1) .ne. ' ') file = temp
        else if ( nfile .eq. 1 ) then
          file = temp
        else if ( nfile .eq. 2 ) then
          call getarg(ifile(2),temp)
          file = temp
        end if
        iloc = index(file,' ')
        file = file(1:(iloc-1))//char(0)

        if (twod) then
          call wdata2d( Ma, Pr, Re, time, q(1,1,1,1), q(1,1,1,2), &
                        q(1,1,1,3), q(1,1,1,5), nx, ny, nx, file )
        else
          call wdata  ( Ma, Pr, Re, time, q(1,1,1,1), q(1,1,1,2), &
                        q(1,1,1,3), q(1,1,1,4), q(1,1,1,5), &
                        nx, ny, nz, nx, ny, file )
        end if

        end if

        call exit(0)

 1000   write(*,"('>> Error in file I/O')")
        call exit(1)

        end program conv


text:
      subroutine sub(num, nums1, nums2)
      dimension nums1(6)
      dimension nums2(2,3)
      write(6, '(i1)') num
      do i=1,6
        write(6, '(i2)') nums1(i)
      enddo
      do i=1,2
        do j=1,3
          write(6, '(i2)') nums2(i,j)
        enddo
      enddo
      end

      program prog
      dimension nums(6)
      nums(1) = 3
      do i=2,6
        nums(i) = nums(i-1) + i
      enddo
      call sub(nums, nums, nums)
      end