text:
MODULE sweep_module

  use parallel_module
  use electron_module, only: Nfixed, Ndspin, Nspin, Nband, Nelectron
  use bz_module, only: weight_bz, Nbzsm
  use wf_module
  use cg_module, only: conjugate_gradient
  use array_bound_module, only: ML_0,ML_1,MBZ_0,MBZ_1,MSP_0,MSP_1,MB_0,MB_1
  use gram_schmidt_module
  use io_module
  use total_energy_module, only: calc_with_rhoIN_total_energy
  use fermi_module
  use subspace_diag_module
  use esp_gather_module
  use watch_module
  use hamiltonian_module
  use xc_hybrid_module, only: control_xc_hybrid, get_flag_xc_hybrid
  use io_tools_module
  use eigenvalues_module

  implicit none

  PRIVATE
  PUBLIC :: calc_sweep, init_sweep, read_sweep

  integer :: iconv_check=1
  real(8) :: Echk, Echk0
  real(8) :: tol_Echk=1.d-12
  real(8) :: tol_esp=1.d-7
  real(8) :: max_esperr
  integer :: mb_ref
  integer :: Nsweep

CONTAINS

  SUBROUTINE read_sweep
    implicit none
    call IOTools_readIntegerKeyword( "NSWEEP", Nsweep )
  END SUBROUTINE read_sweep

  SUBROUTINE init_sweep( iconv_check_in, mb_ref_in, tol_in )
    implicit none
    integer,intent(IN) :: iconv_check_in, mb_ref_in
    real(8),intent(IN) :: tol_in
    iconv_check = iconv_check_in
    mb_ref = mb_ref_in
    select case( iconv_check )
    case( 1 )
       tol_Echk = tol_in
    case( 2 )
       tol_esp = tol_in
    case( 3 )
       tol_esp = tol_in
    end select
  END SUBROUTINE init_sweep

  SUBROUTINE calc_sweep( disp_switch, ierr_out, Diter_in, outer_loop_info )
    implicit none
    logical,intent(IN)  :: disp_switch
    integer,intent(OUT) :: ierr_out
    integer,optional,intent(IN)  :: Diter_in
    character(*),optional,intent(IN) :: outer_loop_info
    integer :: iter,s,k,n,m,iflag_hybrid,ierr,Diter
    logical :: flag_exit, flag_conv
    logical :: flag_end, flag_end1, flag_end2
    character(40) :: chr_iter
    character(22) :: add_info
    type(time) :: etime
    type(eigv) :: eval
    logical,external :: exit_program

    Diter = 0
    if ( present(Diter_in) ) then
       Diter = Diter_in
    else
       Diter = Nsweep
    end if

    if ( Diter <= 0 ) return

    call write_border( 0, "" )
    call write_border( 0, " SWEEP START ----------" )

    flag_end  = .false.
    flag_exit = .false.
    flag_conv = .false.
    Echk      = 0.0d0
    ierr_out  = 0
    add_info  = "" ; if ( present(outer_loop_info) ) add_info=outer_loop_info

    allocate( esp0(Nband,Nbzsm,Nspin) ) ; esp0=0.0d0

    call get_flag_xc_hybrid( iflag_hybrid )

    select case( iflag_hybrid )
    case(0,1,3)

       if ( iflag_hunk >= 1 ) then
          do s=MSP_0,MSP_1
          do k=MBZ_0,MBZ_1
             do m=MB_0,MB_1,MB_d
                n=min(m+MB_d-1,MB_1)
                call hamiltonian &
                     (k,s,unk(ML_0,m,k,s),hunk(ML_0,m,k,s),ML_0,ML_1,m,n)
             end do
          end do
          end do
       end if

    case( 2 )

       if ( iflag_hunk >= 1 ) then
          call control_xc_hybrid(0)
          allocate( workwf(ML_0:ML_1,MB_d) ) ; workwf=0.0d0
          do s=MSP_0,MSP_1
          do k=MBZ_0,MBZ_1
             do m=MB_0,MB_1,MB_d
                n=min(m+MB_d-1,MB_1)
                workwf(:,1:n-m+1)=hunk(:,m:n,k,s)
                call hamiltonian &
                     (k,s,unk(ML_0,m,k,s),hunk(ML_0,m,k,s),ML_0,ML_1,m,n)
                hunk(:,m:n,k,s)=hunk(:,m:n,k,s)+workwf(:,1:n-m+1)
             end do ! m
          end do ! k
          end do ! s
          deallocate( workwf )
       end if

       call control_xc_hybrid(1)

    end select

    call calc_with_rhoIN_total_energy( Echk )

    do iter=1,Diter

       write(chr_iter,'(" sweep_iter=",i4,1x,a)') iter, add_info
       call write_border( 0, chr_iter(1:len_trim(chr_iter)) )

       call init_time_watch( etime )

       Echk0=Echk
       esp0 =esp
       do s=MSP_0,MSP_1
       do k=MBZ_0,MBZ_1

          call conjugate_gradient( ML_0,ML_1, Nband, k,s &
                                 ,unk(ML_0,1,k,s), esp(1,k,s), res(1,k,s) )

          call gram_schmidt(1,Nband,k,s)

          call subspace_diag(k,s)

       end do
       end do

       call esp_gather(Nband,Nbzsm,Nspin,esp)

#ifdef _DRSDFT_
       call mpi_bcast( unk, size(unk), MPI_REAL8, 0, comm_fkmb, ierr )
#else
       call mpi_bcast( unk, size(unk), MPI_COMPLEX16, 0, comm_fkmb, ierr )
#endif
       call mpi_bcast( esp, size(esp), MPI_REAL8, 0, comm_fkmb, ierr )

       call calc_fermi(iter,Nfixed,Nband,Nbzsm,Nspin,Nelectron,Ndspin &
                      ,esp,weight_bz,occ,disp_switch)

       call calc_with_rhoIN_total_energy( Echk )

       call conv_check( iter, res, flag_conv )
       call global_watch( .false., flag_end1 )
       flag_end2 = exit_program()
       flag_end  = ( flag_end1 .or. flag_end2 )
       flag_exit = (flag_end.or.flag_conv.or.(iter==Diter))

       if ( disp_switch ) call write_info_sweep
       call construct_eigenvalues( Nband, Nbzsm, Nspin, esp, eval )
       if ( myrank == 0 ) call write_eigenvalues( eval )

       call calc_time_watch( etime )
       if ( disp_switch ) then
          write(*,*)
          write(*,'(1x,"time(sweep)=",f10.3,"(rank0)",f10.3,"(min)" &
               ,f10.3,"(max)")') etime%t0, etime%tmin, etime%tmax
          write(*,*)
       end if

       call write_data( disp_switch, flag_exit )

       if ( flag_exit ) exit

    end do ! iter

    if ( myrank == 0 ) call write_info_esp_wf( 2 )

    deallocate( esp0 )

    ierr_out = iter

    if ( flag_end1 ) then
       ierr_out = -1
       if ( myrank == 0 ) write(*,*) "flag_end=",flag_end
       return
    end if

    if ( flag_end2 ) then
       ierr_out = -3
       if ( myrank == 0 ) write(*,*) "flag_end=",flag_end
       return
    end if

    if ( iter > Diter ) then
       ierr_out = -2
       if ( myrank == 0 ) write(*,*) "sweep not converged"
    end if

    call gather_wf

    call write_border( 0, " SWEEP END ----------" )
    call write_border( 0, "" )

  END SUBROUTINE calc_sweep

  SUBROUTINE conv_check( iter, res, flag_conv )
    implicit none
    integer,intent(IN)  :: iter
    real(8),intent(INOUT) :: res(:,:,:)
    logical,intent(OUT) :: flag_conv
    call write_border( 1, " conv_check(start)" )
    select case( iconv_check )
    case( 1 )
         call conv_check_1( iter, flag_conv )
    case( 2 )
         call conv_check_2( flag_conv )
    case( 3 )
         call conv_check_3( res, flag_conv )
    end select
    call write_border( 1, " conv_check(end)" )
  END SUBROUTINE conv_check

  SUBROUTINE conv_check_1( iter, flag_conv )
    implicit none
    integer,intent(IN)  :: iter
    logical,intent(OUT) :: flag_conv
    flag_conv=.false.
    if ( abs(Echk-Echk0) < tol_Echk ) flag_conv=.true.
  END SUBROUTINE conv_check_1

  SUBROUTINE conv_check_2( flag_conv )
    implicit none
    logical,intent(OUT) :: flag_conv
    integer :: ierr
    real(8) :: err0
    max_esperr = maxval( abs(  esp(1:mb_ref,MBZ_0,MSP_0:MSP_1) &
                             -esp0(1:mb_ref,MBZ_0,MSP_0:MSP_1) ) )
    call mpi_allreduce(max_esperr,err0,1,MPI_REAL8,MPI_MAX,comm_spin,ierr)
    call mpi_allreduce(err0,max_esperr,1,MPI_REAL8,MPI_MAX,comm_bzsm,ierr)
    flag_conv = .false.
    if ( max_esperr < tol_esp ) flag_conv = .true.
  END SUBROUTINE conv_check_2

  SUBROUTINE conv_check_3( res, flag_conv )
    implicit none
    real(8),intent(OUT) :: res(:,:,:)
    logical,intent(OUT) :: flag_conv
    integer :: ierr
    real(8) :: err0
    res=0.0d0
    res=abs(esp-esp0)
    max_esperr = maxval( res(1:mb_ref,:,:) )
    flag_conv = .false.
    if ( max_esperr < tol_esp ) flag_conv = .true.
    where( res < tol_esp ) res=-1.0d0
  END SUBROUTINE conv_check_3

  SUBROUTINE write_info_sweep
    implicit none
    integer :: s,k,n
    call write_border( 1, " write_info_sweep(start)" )
!    write(*,'(a4,a6,a20,2a13,1x)') &
!         "k","n","esp(n,k,s)","esp_err  ","occ(n,k,s)  "
!    do k=1,Nbzsm
!    do n=max(1,nint(Nelectron/2)-5),min(nint(Nelectron/2)+5,Nband)
!       write(*,'(i4,i6,2(f20.15,2g13.5,1x))') k,n &
!            ,(esp(n,k,s),esp(n,k,s)-esp0(n,k,s),occ(n,k,s),s=1,Nspin)
!    end do
!    end do
    if ( iconv_check == 1 ) then
       write(*,'(/,1x,"Echk,dif/tol =",g18.10,2x,g12.5," /",es12.5)') &
            Echk, Echk-Echk0, tol_Echk
    else
       write(*,'(/,1x,"max_esperr/tol, mb_ref =",g12.5," /",es12.5,i7)') &
            max_esperr,tol_esp,mb_ref
    end if
    call write_esp_wf
!    write(*,*) "sum(occ)=",(sum(occ(:,:,s)),s=1,Nspin)
    call write_border( 1, " write_info_sweep(end)" )
  END SUBROUTINE write_info_sweep

END MODULE sweep_module

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
      MODULE GWFEVTMODULE
        INTEGER, SAVE, POINTER                 ::NEVTOP,IEVTCB
        INTEGER, SAVE, POINTER                 ::NPEVT,IEVTPF
        REAL,    SAVE,   DIMENSION(:,:),  POINTER      ::EVTR
        REAL,    SAVE,   DIMENSION(:,:),  POINTER      ::EXDP
        REAL,    SAVE,   DIMENSION(:,:),  POINTER      ::SURF
        INTEGER, SAVE,   DIMENSION(:,:),  POINTER      ::IEVT
      TYPE GWFEVTTYPE
        INTEGER,  POINTER                 ::NEVTOP,IEVTCB
        INTEGER,  POINTER                 ::NPEVT,IEVTPF
        REAL,       DIMENSION(:,:),  POINTER      ::EVTR
        REAL,       DIMENSION(:,:),  POINTER      ::EXDP
        REAL,       DIMENSION(:,:),  POINTER      ::SURF
        INTEGER,    DIMENSION(:,:),  POINTER      ::IEVT
      END TYPE
      TYPE(GWFEVTTYPE), SAVE ::GWFEVTDAT(10)
      END MODULE GWFEVTMODULE

      SUBROUTINE GWF2EVT7AR(IN,IGRID)
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR EVAPOTRANSPIRATION
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:IOUT,NCOL,NROW,IFREFM
      USE GWFEVTMODULE,ONLY:NEVTOP,IEVTCB,NPEVT,IEVTPF,EVTR,EXDP,SURF,
     1                      IEVT
C
      CHARACTER*200 LINE
      CHARACTER*4 PTYP
C     ------------------------------------------------------------------
C
C1------ALLOCATE SCALAR VARIABLES.
      ALLOCATE(NEVTOP,IEVTCB)
      ALLOCATE(NPEVT,IEVTPF)
C
C2------IDENTIFY PACKAGE.
      IEVTPF=0
      WRITE(IOUT,1)IN
    1 FORMAT(1X,/1X,'EVT -- EVAPOTRANSPIRATION PACKAGE, VERSION 7,',
     1     ' 5/2/2005',/,9X,'INPUT READ FROM UNIT ',I4)
C
C3------READ ET OPTION (NEVTOP) AND UNIT OR FLAG FOR CELL-BY-CELL FLOW
C3------TERMS (IEVTCB).
      CALL URDCOM(IN,IOUT,LINE)
      CALL UPARARRAL(IN,IOUT,LINE,NPEVT)
      IF(IFREFM.EQ.0) THEN
         READ(LINE,'(2I10)') NEVTOP,IEVTCB
      ELSE
         LLOC=1
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NEVTOP,R,IOUT,IN)
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IEVTCB,R,IOUT,IN)
      END IF
C
C4------CHECK TO SEE THAT ET OPTION IS LEGAL.
      IF(NEVTOP.GE.1.AND.NEVTOP.LE.3)GO TO 200
C
C4A-----OPTION IS ILLEGAL -- PRINT A MESSAGE & ABORT SIMULATION.
      WRITE(IOUT,8) NEVTOP
    8 FORMAT(1X,'ILLEGAL ET OPTION CODE (NEVTOP = ',I5,
     1       ') -- SIMULATION ABORTING')
      CALL USTOP(' ')
C
C5------OPTION IS LEGAL -- PRINT THE OPTION CODE.
  200 IF(NEVTOP.EQ.1) WRITE(IOUT,201)
  201 FORMAT(1X,'OPTION 1 -- EVAPOTRANSPIRATION FROM TOP LAYER')
      IF(NEVTOP.EQ.2) WRITE(IOUT,202)
  202 FORMAT(1X,'OPTION 2 -- EVAPOTRANSPIRATION FROM ONE SPECIFIED',
     1   ' NODE IN EACH VERTICAL COLUMN')
      IF(NEVTOP.EQ.3) WRITE(IOUT,203)
  203 FORMAT(1X,'OPTION 3 -- EVAPOTRANSPIRATION FROM HIGHEST ACTIVE',
     1   ' NODE IN EACH VERTICAL COLUMN')
C
C6------IF CELL-BY-CELL FLOWS ARE TO BE SAVED, THEN PRINT UNIT NUMBER.
      IF(IEVTCB.GT.0) WRITE(IOUT,204) IEVTCB
  204 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE SAVED ON UNIT ',I4)
C
C7------ALLOCATE SPACE FOR THE ARRAYS EVTR, EXDP, SURF, AND IEVT.
      ALLOCATE (EVTR(NCOL,NROW))
      ALLOCATE (EXDP(NCOL,NROW))
      ALLOCATE (SURF(NCOL,NROW))
      ALLOCATE (IEVT(NCOL,NROW))
C
C8------READ NAMED PARAMETERS
      WRITE(IOUT,5) NPEVT
    5 FORMAT(1X,//1X,I5,' Evapotranspiration parameters')
      IF(NPEVT.GT.0) THEN
         DO 20 K=1,NPEVT
         CALL UPARARRRP(IN,IOUT,N,0,PTYP,1,1,0)
         IF(PTYP.NE.'EVT') THEN
            WRITE(IOUT,7)
    7       FORMAT(1X,'Parameter type must be EVT')
            CALL USTOP(' ')
         END IF
   20    CONTINUE
      END IF
C
C9------RETURN
      CALL SGWF2EVT7PSV(IGRID)
      RETURN
      END
      SUBROUTINE GWF2EVT7RP(IN,IGRID)
C     ******************************************************************
C     READ EVAPOTRANSPIRATION DATA
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,DELR,DELC,IFREFM
      USE GWFEVTMODULE,ONLY:NEVTOP,NPEVT,IEVTPF,EVTR,EXDP,SURF,IEVT
C
      CHARACTER*24 ANAME(4)
C
      DATA ANAME(1) /'          ET LAYER INDEX'/
      DATA ANAME(2) /'              ET SURFACE'/
      DATA ANAME(3) /' EVAPOTRANSPIRATION RATE'/
      DATA ANAME(4) /'        EXTINCTION DEPTH'/
C     ------------------------------------------------------------------
C
C1------SET POINTERS FOR THE CURRENT GRID.
      CALL SGWF2EVT7PNT(IGRID)
C
C2------READ FLAGS SHOWING WHETHER DATA IS TO BE REUSED.
      IF(NEVTOP.EQ.2) THEN
         IF(IFREFM.EQ.0) THEN
            READ(IN,'(4I10)') INSURF,INEVTR,INEXDP,INIEVT
         ELSE
            READ(IN,*) INSURF,INEVTR,INEXDP,INIEVT
         END IF
      ELSE
         IF(IFREFM.EQ.0) THEN
            READ(IN,'(3I10)') INSURF,INEVTR,INEXDP
         ELSE
            READ(IN,*) INSURF,INEVTR,INEXDP
         END IF
      END IF
C
C3------TEST INSURF TO SEE WHERE SURFACE ELEVATION COMES FROM.
      IF(INSURF.LT.0) THEN
C
C3A------INSURF<0, SO REUSE SURFACE ARRAY FROM LAST STRESS PERIOD.
        WRITE(IOUT,3)
    3   FORMAT(1X,/1X,'REUSING SURF FROM LAST STRESS PERIOD')
      ELSE
C
C3B------INSURF=>0, SO CALL MODULE U2DREL TO READ SURFACE.
        CALL U2DREL(SURF,ANAME(2),NROW,NCOL,0,IN,IOUT)
      END IF
C
C4------TEST INEVTR TO SEE WHERE MAX ET RATE (EVTR) COMES FROM.
      IF(INEVTR.LT.0) THEN
C
C4A-----INEVTR<0, SO REUSE EVTR FROM LAST STRESS PERIOD.
        WRITE(IOUT,4)
    4   FORMAT(1X,/1X,'REUSING EVTR FROM LAST STRESS PERIOD')
      ELSE
C
C4B-----INEVTR=>0, SO READ MAX ET RATE.
        IF(NPEVT.EQ.0) THEN
C
C4B1----THERE ARE NO PARAMETERS,SO READ EVTR USING U2DREL
          CALL U2DREL(EVTR,ANAME(3),NROW,NCOL,0,IN,IOUT)
        ELSE
C
C4B2----DEFINE EVTR USING PARAMETERS. INEVTR IS THE NUMBER OF
C4B2----PARAMETERS TO USE THIS STRESS PERIOD.
          CALL PRESET('EVT')
          WRITE(IOUT,33)
   33     FORMAT(1X,///1X,
     1      'EVTR array defined by the following parameters:')
          IF (INEVTR.EQ.0) THEN
            WRITE(IOUT,34)
   34       FORMAT(' ERROR: When parameters are defined for the EVT',
     &     ' Package, at least one parameter',/,' must be specified',
     &     ' each stress period -- STOP EXECUTION (GWF2EVT7RPSS)')
            CALL USTOP(' ')
          END IF
          CALL UPARARRSUB2(EVTR,NCOL,NROW,0,INEVTR,IN,IOUT,'EVT',
     1        ANAME(3),'EVT',IEVTPF)
        END IF
C
C5------MULTIPLY MAX ET RATE BY CELL AREA TO GET VOLUMETRIC RATE
        DO 40 IR=1,NROW
        DO 40 IC=1,NCOL
        EVTR(IC,IR)=EVTR(IC,IR)*DELR(IC)*DELC(IR)
   40   CONTINUE
      END IF
C
C6------TEST INEXDP TO SEE WHERE EXTINCTION DEPTH COMES FROM
      IF(INEXDP.LT.0) THEN
C
C6A------IF INEXDP<0 REUSE EXTINCTION DEPTH FROM LAST STRESS PERIOD
        WRITE(IOUT,5)
    5   FORMAT(1X,/1X,'REUSING EXDP FROM LAST STRESS PERIOD')
      ELSE
C
C6B------IF INEXDP=>0 CALL MODULE U2DREL TO READ EXTINCTION DEPTH
        CALL U2DREL(EXDP,ANAME(4),NROW,NCOL,0,IN,IOUT)
      END IF
C
C7------IF OPTION(NEVTOP) IS 2 THEN WE NEED AN INDICATOR ARRAY.  TEST
C7------INIEVT TO SEE HOW TO DEFINE IEVT.
      IF(NEVTOP.EQ.2) THEN
        IF(INIEVT.LT.0) THEN
C
C7A------IF INIEVT<0 THEN REUSE LAYER INDICATOR ARRAY.
          WRITE(IOUT,2)
    2     FORMAT(1X,/1X,'REUSING IEVT FROM LAST STRESS PERIOD')
        ELSE
C
C7B------IF INIEVT=>0 THEN CALL MODULE U2DINT TO READ INDICATOR ARRAY.
   49     CALL U2DINT(IEVT,ANAME(1),NROW,NCOL,0,IN,IOUT)
          DO 57 IR=1,NROW
          DO 57 IC=1,NCOL
          IF(IEVT(IC,IR).LT.1 .OR. IEVT(IC,IR).GT.NLAY) THEN
            WRITE(IOUT,56) IC,IR,IEVT(IC,IR)
   56       FORMAT(1X,/1X,'INVALID LAYER NUMBER IN IEVT FOR COLUMN',I4,
     1           '  ROW',I4,'  :',I4)
            CALL USTOP(' ')
          END IF
   57     CONTINUE
        END IF
      END IF
C
C8------RETURN
      RETURN
      END
      SUBROUTINE GWF2EVT7FM(IGRID)
C     ******************************************************************
C     ADD EVAPOTRANSPIRATION TO RHS AND HCOF
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND,HNEW,RHS,HCOF,BOTM
      USE GWFEVTMODULE,ONLY:NEVTOP,EVTR,EXDP,SURF,IEVT
C
      DOUBLE PRECISION HH,SS,XX,DD
C     ------------------------------------------------------------------
C
C1------SET POINTERS FOR THE CURRENT GRID.
      CALL SGWF2EVT7PNT(IGRID)
C
C2------PROCESS EACH HORIZONTAL CELL LOCATION
      DO 10 IR=1,NROW
      DO 10 IC=1,NCOL
C
C3------SET THE LAYER INDEX -- FOR OPTION 1, THE LAYER IS 1;
C3------FOR OPTION 2, THE LAYER IS SPECIFIED IN IEVT.
      IF(NEVTOP.EQ.1) THEN
         IL=1
      ELSE IF(NEVTOP.EQ.2) THEN
         IL=IEVT(IC,IR)
         IF(IL.EQ.0) GO TO 10
      ELSE
C
C4------FOR OPTION 3, FIND UPPERMOST ACTIVE CELL.
         DO 3 IL=1,NLAY
         IF(IBOUND(IC,IR,IL).NE.0) GO TO 4
    3    CONTINUE
         IL=1
      END IF
C
C5------IF THE CELL IS NOT VARIABLE HEAD, IGNORE IT.  IF CELL IS
C5------VARIABLE HEAD, GET DATA NEEDED TO COMPUTE FLOW TERMS.
    4 IF(IBOUND(IC,IR,IL).LE.0)GO TO 10
      C=EVTR(IC,IR)
      S=SURF(IC,IR)
      SS=S
      HH=HNEW(IC,IR,IL)
C
C6------IF AQUIFER HEAD IS GREATER THAN OR EQUAL TO SURF, ET IS CONSTANT
      IF(HH.LT.SS) GO TO 5
C
C6A-----HEAD IS GREATER THAN OR EQUAL TO SURF.  ADD EVTR TO RHS
      RHS(IC,IR,IL)=RHS(IC,IR,IL) + C
      GO TO 10
C
C7------IF DEPTH TO WATER>=EXTINCTION DEPTH, THEN ET IS 0.
    5 DD=SS-HH
      X=EXDP(IC,IR)
      XX=X
      IF(DD.GE.XX)GO TO 10
C
C8------LINEAR RANGE. ADD ET TERMS TO BOTH RHS AND HCOF.
      RHS(IC,IR,IL)=RHS(IC,IR,IL)+C-C*S/X
      HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-C/X
   10 CONTINUE
C
C9------RETURN
      RETURN
      END
      SUBROUTINE GWF2EVT7BD(KSTP,KPER,IGRID)
C     ******************************************************************
C     CALCULATE VOLUMETRIC BUDGET FOR EVAPOTRANSPIRATION
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IBOUND,HNEW,BUFF
      USE GWFBASMODULE,ONLY:MSUM,VBVL,VBNM,ICBCFL,DELT,PERTIM,TOTIM
      USE GWFEVTMODULE,ONLY:NEVTOP,IEVTCB,EVTR,EXDP,SURF,IEVT
C
      DOUBLE PRECISION RATOUT,QQ,HH,SS,DD,XX,HHCOF,RRHS
      CHARACTER*16 TEXT
      DATA TEXT /'              ET'/
C     ------------------------------------------------------------------
C
C1------SET POINTERS FOR THE CURRENT GRID.
      CALL SGWF2EVT7PNT(IGRID)
C
C2------CLEAR THE RATE ACCUMULATOR.
      ZERO=0.
      RATOUT=ZERO
C
C3------CLEAR THE BUFFER & SET CELL-BY-CELL BUDGET SAVE FLAG (IBD).
      DO 2 IL=1,NLAY
      DO 2 IR=1,NROW
      DO 2 IC=1,NCOL
      BUFF(IC,IR,IL)=ZERO
    2 CONTINUE
      IBD=0
      IF(IEVTCB.GT.0) IBD=ICBCFL
C
C4------PROCESS EACH HORIZONTAL CELL LOCATION.
      DO 10 IR=1,NROW
      DO 10 IC=1,NCOL
C
C5------SET THE LAYER INDEX -- FOR OPTION 1, THE LAYER IS 1;
C5------FOR OPTION 2, THE LAYER IS SPECIFIED IN IEVT.
      IF(NEVTOP.EQ.1) THEN
         IL=1
      ELSE IF(NEVTOP.EQ.2) THEN
         IL=IEVT(IC,IR)
         IF(IL.EQ.0) GO TO 10
      ELSE
C
C6------FOR OPTION 3, FIND UPPERMOST NON-DRY CELL.
         DO 3 IL=1,NLAY
         IF(IBOUND(IC,IR,IL).NE.0) GO TO 4
    3    CONTINUE
         IL=1
    4    IEVT(IC,IR)=IL
      END IF
C
C7------IF CELL IS EXTERNAL THEN IGNORE IT.
      IF(IBOUND(IC,IR,IL).LE.0)GO TO 10
      C=EVTR(IC,IR)
      S=SURF(IC,IR)
      SS=S
      HH=HNEW(IC,IR,IL)
C
C8------IF AQUIFER HEAD => SURF,SET Q=MAX ET RATE.
      IF(HH.LT.SS) GO TO 7
      QQ=-C
      GO TO 9
C
C9------IF DEPTH=>EXTINCTION DEPTH, ET IS 0.
    7 X=EXDP(IC,IR)
      XX=X
      DD=SS-HH
      IF(DD.GE.XX)GO TO 10
C
C10-----LINEAR RANGE. Q= -HNEW*EVTR/EXDP -EVTR + EVTR*SURF/EXDP.
      HHCOF=-C/X
      RRHS=(C*S/X)-C
      QQ=HH*HHCOF+RRHS
C
C11-----ACCUMULATE TOTAL FLOW RATE.
    9 Q=QQ
      RATOUT=RATOUT-QQ
C
C12-----ADD Q TO BUFFER.
      BUFF(IC,IR,IL)=Q
   10 CONTINUE
C
C13-----IF CELL-BY-CELL FLOW TO BE SAVED, CALL APPROPRIATE UTILITY
C13-----MODULE SAVE THEM.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IEVTCB,BUFF,NCOL,NROW,
     1                          NLAY,IOUT)
      IF(IBD.EQ.2) CALL UBDSV3(KSTP,KPER,TEXT,IEVTCB,BUFF,IEVT,NEVTOP,
     1                   NCOL,NROW,NLAY,IOUT,DELT,PERTIM,TOTIM,IBOUND)
C
C14-----MOVE TOTAL ET RATE INTO VBVL FOR PRINTING BY BAS1OT.
      ROUT=RATOUT
      VBVL(3,MSUM)=ZERO
      VBVL(4,MSUM)=ROUT
C
C15-----ADD ET(ET_RATE TIMES STEP LENGTH) TO VBVL.
      VBVL(2,MSUM)=VBVL(2,MSUM)+ROUT*DELT
C
C16-----MOVE BUDGET TERM LABELS TO VBNM FOR PRINT BY MODULE BAS1OT.
      VBNM(MSUM)=TEXT
C
C17-----INCREMENT BUDGET TERM COUNTER.
      MSUM=MSUM+1
C
C18-----RETURN.
      RETURN
      END
      SUBROUTINE GWF2EVT7DA(IGRID)
C  Deallocate EVT MEMORY
      USE GWFEVTMODULE
C
        DEALLOCATE(GWFEVTDAT(IGRID)%NEVTOP)
        DEALLOCATE(GWFEVTDAT(IGRID)%IEVTCB)
        DEALLOCATE(GWFEVTDAT(IGRID)%NPEVT)
        DEALLOCATE(GWFEVTDAT(IGRID)%IEVTPF)
        DEALLOCATE(GWFEVTDAT(IGRID)%EVTR)
        DEALLOCATE(GWFEVTDAT(IGRID)%EXDP)
        DEALLOCATE(GWFEVTDAT(IGRID)%SURF)
        DEALLOCATE(GWFEVTDAT(IGRID)%IEVT)
C
      RETURN
      END
      SUBROUTINE SGWF2EVT7PNT(IGRID)
C  Set pointers to EVT data for grid.
      USE GWFEVTMODULE
C
        NEVTOP=>GWFEVTDAT(IGRID)%NEVTOP
        IEVTCB=>GWFEVTDAT(IGRID)%IEVTCB
        NPEVT=>GWFEVTDAT(IGRID)%NPEVT
        IEVTPF=>GWFEVTDAT(IGRID)%IEVTPF
        EVTR=>GWFEVTDAT(IGRID)%EVTR
        EXDP=>GWFEVTDAT(IGRID)%EXDP
        SURF=>GWFEVTDAT(IGRID)%SURF
        IEVT=>GWFEVTDAT(IGRID)%IEVT
C
      RETURN
      END
      SUBROUTINE SGWF2EVT7PSV(IGRID)
C  Save pointers to EVT data for grid.
      USE GWFEVTMODULE
C
        GWFEVTDAT(IGRID)%NEVTOP=>NEVTOP
        GWFEVTDAT(IGRID)%IEVTCB=>IEVTCB
        GWFEVTDAT(IGRID)%NPEVT=>NPEVT
        GWFEVTDAT(IGRID)%IEVTPF=>IEVTPF
        GWFEVTDAT(IGRID)%EVTR=>EVTR
        GWFEVTDAT(IGRID)%EXDP=>EXDP
        GWFEVTDAT(IGRID)%SURF=>SURF
        GWFEVTDAT(IGRID)%IEVT=>IEVT
C
      RETURN
      END

ext: f

max_stars_repo_name: cneyens/MODFLOW-MAXDRN

max_issues_repo_name: cneyens/MODFLOW-MAXDRN

max_forks_repo_name: cneyens/MODFLOW-MAXDRN


text:
*> \brief \b SNRM2
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       REAL FUNCTION SNRM2(N,X,INCX)
*
*       .. Scalar Arguments ..
*       INTEGER INCX,N
*       ..
*       .. Array Arguments ..
*       REAL X(*)
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SNRM2 returns the euclidean norm of a vector via the function
*> name, so that
*>
*>    SNRM2 := sqrt( x'*x ).
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup single_blas_level1
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  -- This version written on 25-October-1982.
*>     Modified on 14-October-1993 to inline the call to SLASSQ.
*>     Sven Hammarling, Nag Ltd.
*> \endverbatim
*>
*  =====================================================================
      REAL FUNCTION SNRM2(N,X,INCX)
*
*  -- Reference BLAS level1 routine (version 3.4.0) --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER INCX,N
*     ..
*     .. Array Arguments ..
      REAL X(*)
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL ONE,ZERO
      PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
*     ..
*     .. Local Scalars ..
      REAL ABSXI,NORM,SCALE,SSQ
      INTEGER IX
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SQRT
*     ..
      IF (N.LT.1 .OR. INCX.LT.1) THEN
          NORM = ZERO
      ELSE IF (N.EQ.1) THEN
          NORM = ABS(X(1))
      ELSE
          SCALE = ZERO
          SSQ = ONE
*        The following loop is equivalent to this call to the LAPACK
*        auxiliary routine:
*        CALL SLASSQ( N, X, INCX, SCALE, SSQ )
*
          DO 10 IX = 1,1 + (N-1)*INCX,INCX
              IF (X(IX).NE.ZERO) THEN
                  ABSXI = ABS(X(IX))
                  IF (SCALE.LT.ABSXI) THEN
                      SSQ = ONE + SSQ* (SCALE/ABSXI)**2
                      SCALE = ABSXI
                  ELSE
                      SSQ = SSQ + (ABSXI/SCALE)**2
                  END IF
              END IF
   10     CONTINUE
          NORM = SCALE*SQRT(SSQ)
      END IF
*
      SNRM2 = NORM
      RETURN
*
*     End of SNRM2.
*
      END

max_stars_repo_name: haampie/libflame

max_issues_repo_name: haampie/libflame

max_forks_repo_name: haampie/libflame


text:
 General information
A wiki:WikiPedia:Didgeridoo didjeridu is a hollow pipe typically 46 feet long and around 2 inches in diameter, though a wide range of sizes exist. Traditional didjeridus are most notably made from eucalyptus limbs partially hollowed out by termites. There is also documentation of traditional didjeridus made from palm, and it is believed that the first ones were made from bamboo. Usually a mouthpiece of wax or hardened gum is on one end and the player blows into the instrument in a trumpet like fashion. The didjeridu is part of the 423.121.11 group in the wiki:WikiPedia:Aerophone aerophone class from the wiki:WikiPedia:HornbostelSachs HornbostelSachs system of musical classification. When applying a western orchestral based classification, the didjeridu is classified as a wiki:WikiPedia:Wind_instrument wind instrument or more specifically a wiki:WikiPedia:Brass_instrument brass instrument. The instrument originated from the indigenous people of Australia, and they call it the yidaki. It is arguably the oldest known musical wind instrument, dating back at least 40,000 years.

Didjeridu has a variety of names/spellings, including:
   didjeridu
   didjeridoo
   didgeridu
   didgeridoo
   didj
   yirdaki
   yidaki
   yedaki

Didjeridus are often associated with newage Hippies hippie type people that like to use them for things like spiritual sound healings.

 Obtaining a didjeridu
 Buying a didjeridu
    Watermelon Music  They typically carry a couple of didjs, usually found near the cash register.
    Whole Earth Festival  Every year theres a booth at WEF that just sells didjeridus.

 Making your own
The didjeridu is a very simple instrument. In essence it is just a hollow tube, and thus many simple, yet functional, didjeridus can be made with minimal supplies and effort. Some say the secret to a good instrument is having the tube taper towards the mouthpiece, that is, not using a pipe which has parallel internal walls.

http://www.rdrop.com/users/mulara/learner.html Instructions for building a plastic or bamboo didjeridu

 Construction Supplies
    PVC pipe from a hardware store such as ACE Hardware ACE
    Bamboo
    Fluorescent light tube guards

 Playing
Many didjeridu players have learned the method known as circular breathing in order to keep the didjeridu buzzing without stopping for a breath. The trick is to maintain air pressure, without blowing, by using your mouth like the airbag of a bagpipe while you inhale through your nose.

 Didjeridu classes
  The Experimental College often offers a http://asucd.ucdavis.edu/experimentalcollege/index.cfm?pageSearch&search_init1&group_id134 Beginning Didgeridoo class.

Users/JabberWokky has appeared on national television playing the didjeridoo as part of prePhantom Menace coverage.  Having stood in line for a week, NBC filmed him and a few others in the http://www.subgenius.com/ SubGenius band http://www.imdb.com/title/tt0077247/ Ralphus and the Didjeri Donts performing http://www.snpp.com/guides/lyrics_list.html Luke be a Jedi Tonight.

20050210 21:21:54 nbsp  The Experimental College has digeridoo classes taught by Min Yang.  I had a great time taking it. Users/BrianneCrabtree

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module Mod_plcd_UPFormulation
   use typre
   use Mod_plcd_BaseElmope
   use Mod_PointerSetter
   use Mod_plcd_ExternalForces
   implicit none
   private
   public SetPointersUPFormulation, SetPointersUPFormulationEndite

   type, extends(PointerSetter) :: SPUPFormulation
contains
      procedure :: SpecificSet => SpecificSetPointersUPFormulation
   end type
   type(SPUPFormulation) :: SetPointersUPFormulation

   type, extends(PointerSetter) :: SPUPFormulationEndite
contains
      procedure :: SpecificSet => SpecificSetPointersUPFormulationEndite
   end type
   type(SPUPFormulationEndite) :: SetPointersUPFormulationEndite

   real(rp), allocatable :: elpre(:)

   real(rp), allocatable :: elresp(:,:),UPResidualProjection(:,:),gradP(:), elMomentumResidual(:,:)
   real(rp) ::  gppre(1)

contains

   !----------------------------------------------------------------------------
   !Setting Pointers

   subroutine SpecificSetPointersUPFormulation(d)
      implicit none
      class(SPUPFormulation) :: d
            !-------------------------------------------------------
            !UPFormulation
            if (a%UseUpFormulation) then
               call ConcatenateProcedures(ProcHook%Initializations,AllocElmope)
               call ConcatenateProcedures(ProcHook%Finalizations,DeallocElmope)
               call ConcatenateProcedures(ProcHook%PreGauss,GathersElmope)
               call ConcatenateProcedures(ProcHook%InGaussElmats,ComputeAndAssemblyGMatrix)
            endif
   end subroutine

   !For actual computations
   !------------------------------------------------------
   subroutine AllocElmope
      call a%Memor%alloc(e%ndime,e%mnode,elMomentumResidual,'elMomentumResidual','UPFormulation')

   end subroutine

   subroutine DeallocElmope
      call a%Memor%dealloc(e%ndime,e%mnode,elMomentumResidual,'elMomentumResidual','UPFormulation')

   end subroutine

   subroutine GathersElmope
      call e%gather(e%ndime,elMomentumResidual,a%UPResidualProjection)

   end subroutine

   subroutine ComputeAndAssemblyGMatrix
      real(rp) :: gelmat(e%ndime,e%pnode,1,e%pnode)

      integer(ip) :: idime,inode,jnode
      real(rp) :: invK,gpMomentumResidual(e%ndime)

      real(rp) ::   tau, G

      !(div v, p)
      forall(idime = 1:e%ndime, inode = 1:e%pnode, jnode = 1:e%pnode)
         gelmat(idime,inode,1,jnode) = e%cartd(idime,inode)*e%shape(jnode,e%igaus)*dvol
      end forall
      !G mat
      GaussElmat(1:e%ndime,1:e%pnode,e%ndime+1:e%ndime+1,1:e%pnode) = GaussElmat(1:e%ndime,1:e%pnode,e%ndime+1:e%ndime+1,1:e%pnode) - gelmat

      !Gt mat
      forall(idime = 1:e%ndime, inode = 1:e%pnode, jnode = 1:e%pnode)
         GaussElmat(e%ndime+1,jnode,idime,inode) = GaussElmat(e%ndime+1,jnode,idime,inode) + gelmat(idime,inode,1,jnode)
      end forall

      !p p mat, without stabilization
      call ElementMatData%GetInverseVolumetricDeformationModulus(e%igaus,invK)

      forall (inode = 1:e%pnode,jnode = 1:e%pnode)
         GaussElmat(e%ndime+1,inode,e%ndime+1,jnode) = GaussElmat(e%ndime+1,inode,e%ndime+1,jnode) + invK*e%shape(inode,e%igaus)*e%shape(jnode,e%igaus)*dvol
      end forall

      !Stabilization contribution
      call ElementMatData%GetSecantCNorm(e%igaus,G)
      call e%elmlen
      tau = a%staco(1)*e%hleng(1)**2/(2*G)

      !qp
      forall (inode = 1:e%pnode,jnode = 1:e%pnode)
         GaussElmat(e%ndime+1,inode,e%ndime+1,jnode) = GaussElmat(e%ndime+1,inode,e%ndime+1,jnode) + tau*dot_product(e%cartd(:,inode),e%cartd(:,jnode))*dvol
      end forall

      !Orthogonal subscales split oss
      !Nothing to be done here because a Newton Raphson scheme is used, it is already in the external forces

   end subroutine

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!111
   !ENDITE
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
   subroutine SpecificSetPointersUPFormulationEndite(d)
      implicit none
      class(SPUPFormulationEndite) :: d
            !-------------------------------------------------------
            !UPFormulation
            if (a%UseUpFormulation) then
               call ConcatenateProcedures(ProcHook%Initializations,Alloc)
               call ConcatenateProcedures(ProcHook%Finalizations,Dealloc)
               call ConcatenateProcedures(ProcHook%InGauss,PressureInterpolatePAndPassItToEMD)
               call ConcatenateProcedures(ProcHook%PreGauss,Gathers)
               call ConcatenateProcedures(ProcHook%InGaussElmats,PressureContributionToInternalForces)
            endif
   end subroutine

   !For actual computations
   !------------------------------------------------------
   subroutine Alloc
      call a%Memor%alloc(e%mnode,elpre,'elpre','UPFormulation')
      call a%Memor%alloc(e%ndime,e%mnode,elMomentumResidual,'elMomentumResidual','UPFormulation')
      call a%Memor%alloc(e%ndime,gradP,'gradP','AllocOSS')

   end subroutine

   subroutine DeAlloc
      implicit none
      integer(ip) :: ipoin,npoin

      !Modify external forces with pressure boundary conditions
      call a%Mesh%GetNpoin(npoin)
      do ipoin = 1,npoin
         if (ipoin == a%nodpr) then
            if (a%UseUPFormulation .and. (a%kfl_confi == 1)) then
               a%ExternalForcesVector(a%ndofbc+1,ipoin) = a%InternalForcesVector(a%ndofbc+1,ipoin)
            endif
         endif
      enddo

      call a%Memor%dealloc(e%mnode,elpre,'elpre','UPFormulation')
      call a%Memor%dealloc(e%ndime,e%mnode,elMomentumResidual,'elMomentumResidual','UPFormulation')
      call a%Memor%dealloc(e%ndime,gradP,'gradP','AllocOSS')
   end subroutine

   subroutine Gathers
      call e%gather(1,elpre,a%Pressure(:,1))
      call e%gather(e%ndime,elMomentumResidual,a%UPResidualProjection)
   end subroutine

   subroutine PressureInterpolatePAndPassItToEMD
      call e%interpg(1,elpre,gppre)

      call ElementMatData%SetPressureForComputeHistoryAndConstitutiveTensor(gppre(1))
   end subroutine

   subroutine PressureContributionToInternalForces
      real(rp) :: gelmat(e%ndime,e%pnode,1,e%pnode), elext(a%ndofn)

      integer(ip) :: idime,inode,jnode
      real(rp) :: invK

      real(rp) :: divu,gpMomentumResidual(e%ndime),DisplacementSubscales(e%ndime)

      real(rp) ::  tau, G

      !Pressure gradient at gauss point
      call e%gradient(1,elpre,gradP)

      !Momentum equation G
      GaussElInternalForces(1:e%ndime,1:e%pnode) = GaussElInternalForces(1:e%ndime,1:e%pnode) - e%cartd(:,1:e%pnode)*gppre(1)*dvol

      !Continuity equation G
      divu = 0.0_rp
      do idime = 1,e%ndime
         divu = divu + gradDisp(idime,idime)
      enddo
      GaussElInternalForces(e%ndime+1,1:e%pnode) = GaussElInternalForces(e%ndime+1,1:e%pnode) + divu*e%shape(1:e%pnode,e%igaus)*dvol

      !p p mat, without stabilization
      call ElementMatData%GetInverseVolumetricDeformationModulus(e%igaus,invK)
      GaussElInternalForces(e%ndime+1,1:e%pnode) = GaussElInternalForces(e%ndime+1,1:e%pnode) + invK*e%shape(1:e%pnode,e%igaus)*gppre(1)*dvol

      !Stabilization contribution
      call ElementMatData%GetSecantCNorm(e%igaus,G)
      call e%elmlen
      tau = a%staco(1)*e%hleng(1)**2/(2*G)

      !Split OSS, velocity subscales
      call e%interpg(e%ndime,elMomentumResidual,gpMomentumResidual)
      !gpMomentumResidual = 0.0_rp

      call GetExternalForces(e,a,elext)

      DisplacementSubscales = tau*(gradP-elext(1:e%ndime)-gpMomentumResidual)
      !qp stabilization
      forall (inode = 1:e%pnode)
         GaussElInternalForces(e%ndime+1,inode) = GaussElInternalForces(e%ndime+1,inode) + dot_product(e%cartd(:,inode),DisplacementSubscales)*dvol
      end forall

      !Store the subscales for postprocessing
      if (a%UPStoreSubscales) then
         a%UPSubscales(ielem)%a(1:e%ndime,e%igaus) = DisplacementSubscales
      endif

   end subroutine

end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
      SUBROUTINE SETLVL (NEWNM,NUMB,OLDNMS,ITEST,IBIT)
C
C     CREATES A NEW SUBSTRUCTURE NEWNM WHERE
C     - NEWNM IS AN INDEPENDENT SUBSTRUCTURE IF NUMB = 0
C     - NEWNM IS REDUCED FROM THE FIRST SUBSTRUCTURE IN THE ARRAY OLDNMS
C     - NEWNM RESULTS FROM COMBINING THE FIRST I SUBSTRUCTURES IN THE
C       ARRAY OLDNMS IF NUMB = I
C
C     THE OUTPUT VARIABLE ITEST TAKES ON ONE OF THE FOLLOWING VALUES
C          4  IF ONE  OR MORE SUBSTRUCTURES IN OLDNMS DO NOT EXIST
C          7  IF NEWNM ALREADY EXISTS
C          8  IF ONE OF THE SUBSTRUCTURES IN OLDNMS HAS ALREADY
C             BEEN USED IN A REDUCTION OR COMBINATION
C          1  OTHERWISE
C
C     IF ITEST IS SET TO 4, NUMB WILL BE SET TO THE NUMBER OF
C     SUBSTRUCTURES IN OLDNMS THAT DO NOT EXIST AND THE FIRST NUMB NAMES
C     IN OLDNMS WILL BE SET TO THE NAMES OF THOSE SUBSTRUCTURES THAT DO
C     NOT EXIST.  BIT IBIT OF THE FIRST MDI WORD IS SET TO INDICATE THE
C     APPROPRIATE TYPE OF SUBSTRUCTURE. IF IBIT IS ZERO NO CHANGE IS
C     MADE TO THE MDI
C
      IMPLICIT INTEGER (A-Z)
      EXTERNAL        LSHIFT,ANDF,ORF,COMPLF
      LOGICAL         DITUP,MDIUP
      DIMENSION       NEWNM(2),OLDNMS(14),IOLD(7),NMSBR(2)
      COMMON /ZZZZZZ/ BUF(1)
      COMMON /SOF   / DIT,DITPBN,DITLBN,DITSIZ,DITNSB,DITBL,
     1                IODUM(8),MDI,MDIPBN,MDILBN,MDIBL,
     2                NXTDUM(15),DITUP,MDIUP
      DATA    IEMPTY/ 4H    /, NMSBR / 4HSETL,4HVL  /
      DATA    LL,CS , HL    /  2,2,2 /
      DATA    IB    / 1     /
C
      CALL CHKOPN (NMSBR(1))
      ITEST = 1
      CALL FDSUB (NEWNM(1),I)
      IF (I  .NE.  -1) GO TO 500
      IF (NUMB .EQ. 0) GO TO 20
C
C     MAKE SURE THAT ALL THE SUBSTRUCTURES IN OLDNMS DO EXIST.
C
      ICOUNT = 0
      DO 10 I = 1,NUMB
      K = 2*(I-1) + 1
      CALL FDSUB (OLDNMS(K),IOLD(I))
      IF (IOLD(I) .GT. 0) GO TO 10
      ICOUNT = ICOUNT + 1
      KK = 2*(ICOUNT-1) + 1
      OLDNMS(KK  ) = OLDNMS(K  )
      OLDNMS(KK+1) = OLDNMS(K+1)
   10 CONTINUE
      IF (ICOUNT .EQ. 0) GO TO 20
      NUMB = ICOUNT
      GO TO 510
   20 CALL CRSUB (NEWNM(1),INEW)
      IF (NUMB .EQ. 0) RETURN
C
C     NEWNM IS NOT A BASIC SUBSTRUCTURE (LEVEL 0).
C     UPDATE NEWNM S DIRECTORY IN THE MDI.
C
      CALL FMDI (INEW,IMDI)
      LLMASK = COMPLF(LSHIFT(1023,20))
      BUF(IMDI+LL) = ORF(ANDF(BUF(IMDI+LL),LLMASK),LSHIFT(IOLD(1),20))
      IF (IBIT .NE. 0) BUF(IMDI+IB) = ORF(BUF(IMDI+IB),LSHIFT(1,IBIT))
      MDIUP = .TRUE.
C
C     UPDATE IN THE MDI THE DIRECTORIES OF THE SUBSTRUCTURES IN OLDNMS.
C
      IF (NUMB .GT. 7) NUMB = 7
      MASKCS = COMPLF(LSHIFT(1023,10))
      DO 50 I = 1,NUMB
      CALL FMDI (IOLD(I),IMDI)
      IF (ANDF(BUF(IMDI+HL),1023) .EQ. 0) GO TO 40
      ICOUNT = I
      GO TO 520
   40 BUF(IMDI+HL) = ORF(BUF(IMDI+HL),INEW)
      MDIUP = .TRUE.
      IF (NUMB .EQ. 1) RETURN
      IF (I .EQ. NUMB) GO TO 130
      BUF(IMDI+CS) = ORF(ANDF(BUF(IMDI+CS),MASKCS),LSHIFT(IOLD(I+1),10))
   50 CONTINUE
  130 BUF(IMDI+CS) = ORF(ANDF(BUF(IMDI+CS),MASKCS),LSHIFT(IOLD(1),10))
      RETURN
C
C     NEWNM ALREADY EXISTS.
C
  500 ITEST = 7
      RETURN
C
C     ONE OR MORE OF THE SUBSTRUCTURES IN OLDNMS DO NOT EXIST.
C
  510 ITEST = 4
      RETURN
C
C     ONE OF THE SUBSTRUCTURES IN OLDNMS HAS ALREADY BEEN USED IN A
C     REDUCTION OR COMBINATION.  REMOVE ALL CHANGES THAT HAVE BEEN MADE.
C
  520 ITEST = 8
      CALL FDIT (INEW,IDIT)
      BUF(IDIT  ) = IEMPTY
      BUF(IDIT+1) = IEMPTY
      DITUP = .TRUE.
      IF (2*INEW .NE. DITSIZ) GO TO 525
      DITSIZ = DITSIZ - 2
  525 DITNSB = DITNSB - 1
      CALL FMDI (INEW,IMDI)
      BUF(IMDI+LL) = ANDF(BUF(IMDI+LL),LLMASK)
      MDIUP  = .TRUE.
      ICOUNT = ICOUNT - 1
      IF (ICOUNT .LT. 1) RETURN
      DO 530 I = 1,ICOUNT
      CALL FMDI (IOLD(I),IMDI)
      BUF(IMDI+HL) = ANDF(BUF(IMDI+HL),COMPLF(1023))
      BUF(IMDI+CS) = ANDF(BUF(IMDI+CS),MASKCS)
      MDIUP = .TRUE.
  530 CONTINUE
      RETURN
      END


text:
      SUBROUTINE ANSWER(A,S,XMEAN,SIGMA,IDX,PHI,L,M,MM,PF,NDX,ADX)
C------- PRINT INTERMEDIATE RESULTS OF REGRESSION ANALYSIS (SWMREG) ---
      CHARACTER*8 ADX
      INTEGER*4 IDX(4)
      REAL*4 A(7,7),S(4,4),XMEAN(4),SIGMA(4),B(4),BSE(4),PF(3)
C-----------------------------------------------------------------------
      DO 410 I=1,MM
      WRITE(8,400) (A(I,J),J=1,MM)
  400 FORMAT(7E18.8)
  410 CONTINUE
      FVE=1.-A(M,M)
      B0=XMEAN(M)
      YSE=77.7
      IF (PHI .GE. 1) YSE=SIGMA(M)*SQRT(ABS(A(M,M)/PHI))
      DO  5 I=1,L
      IF (IDX(I).EQ.0) GO TO  5
      B(I)=A(I,M)* SQRT(ABS(S(M,M)/S(I,I)))
      BSE(I)=YSE* SQRT(ABS(A(I+M,I+M)/S(I,I)))
      B0=B0-B(I)*XMEAN(I)
    5 CONTINUE
      WRITE(8,10) ADX,NDX,FVE,YSE,B0
   10 FORMAT(/,' VARIABLE ', A8, '................',I5
     2,      /,' FRACTION OF VARIATION EXPLAINED..',E18.8
     3,      /,' STANDARD ERROR OF Y..............',E18.8
     4,      /,' CONSTANT IN REGRESSION EQUATION..',E18.8)
      WRITE(8,20)
   20 FORMAT(/,' VARIABLE     COEFFICIENT      STANDARD ERROR'
     1,'     PARTIAL F-VALUE')
      DO 40 I=1,L
      IF (IDX(I).EQ.0) GO TO 40
      WRITE(8,30) I,B(I),BSE(I),PF(I)
   30 FORMAT(I5,3E20.6)
   40 CONTINUE
      RETURN
      END
      SUBROUTINE AZWTOS(DX,DY,NR,WT,KDX,AZ,TEMP,KEY,INS,IEW)
C------- AZIMUTHAL WEIGHTING OF STATIONS BY QUADRANTS ------------------
	include 'H_param.f'
      CHARACTER*1 INS(N_s),IEW(N_s)
      INTEGER*4 KTX(4),KEMP(N_d),KEY(N_d),KDX(N_d)
      REAL*4 TX(4),TXN(4),DX(N_d),DY(N_d),WT(N_d),AZ(N_d),TEMP(N_d)
C-----------------------------------------------------------------------
      J=0
      DO 10 I=1,NR
      IF (WT(I) .EQ. 0.) GO TO 10
      DXI=DX(I)
      DYI=DY(I)
      IF ((DXI.EQ.0.).AND.(DYI.EQ.0.)) GO TO 6
      JI=KDX(I)
      IF (INS(JI) .EQ. 'S') DYI=-DYI
      IF (IEW(JI) .EQ. 'W') DXI=-DXI
      IF (IEW(1).EQ.'E') THEN
      IF (IEW(JI) .EQ. 'W') THEN
      AZ(I)=AMOD(ATAN2(DXI,DYI)*57.29578 - 360., 360.)
      AZ(I)=-AZ(I)
      ELSE
      AZ(I)=AMOD(ATAN2(DXI,DYI)*57.29578 + 360., 360.)
      GOTO 7
      ENDIF
      ENDIF
      AZ(I)=AMOD(ATAN2(DXI,DYI)*57.29578, 360.)
      IF (IEW(JI) .EQ. 'W') THEN
      IF (AZ(I).LT.0) THEN
      AZ(I)=-AZ(I)
      ELSE
      AZ(I)=AMOD(ATAN2(DXI,DYI)*57.29578 - 360., 360.)
      AZ(I)=-AZ(I)
      ENDIF
      ENDIF
      GO TO 7

      GO TO 7
    6 AZ(I)=999.
    7 J=J+1
      TEMP(J)=AZ(I)
   10 CONTINUE
      CALL SORT(TEMP,KEY,J)
      GAP=TEMP(1)+360.-TEMP(J)
      IG=1
      DO 20 I=2,J
      DTEMP=TEMP(I)-TEMP(I-1)
      IF (DTEMP .LE. GAP) GO TO 20
      GAP=DTEMP
      IG=I
   20 CONTINUE
      TX(1)=TEMP(IG)-0.5*GAP
      TX(2)=TX(1)+90.
      TX(3)=TX(1)+180.
      TX(4)=TX(1)+270.
      DO 124 I=1,4
      TXN(I)=0.
      IF (TX(I) .LT. 0.) TX(I)=TX(I)+360.
      IF (TX(I).GT.360.) TX(I)=TX(I)-360.
  124 CONTINUE
      CALL SORT(TX,KTX,4)
      DO 130 I=1,NR
      IF (WT(I) .EQ. 0.) GO TO 130
      IF (AZ(I) .GT. TX(1)) GO TO 126
  125 TXN(1)=TXN(1)+1.
      KEMP(I)=1
      GO TO 130
  126 IF (AZ(I) .GT. TX(2)) GO TO 127
      TXN(2)=TXN(2)+1.
      KEMP(I)=2
      GO TO 130
  127 IF (AZ(I) .GT. TX(3)) GO TO 128
      TXN(3)=TXN(3)+1.
      KEMP(I)=3
      GO TO 130
  128 IF (AZ(I) .GT. TX(4)) GO TO 125
      TXN(4)=TXN(4)+1.
      KEMP(I)=4
  130 CONTINUE
      XN=4
      IF (TXN(1).EQ.0.) XN=XN-1
      IF (TXN(2).EQ.0.) XN=XN-1
      IF (TXN(3).EQ.0.) XN=XN-1
      IF (TXN(4).EQ.0.) XN=XN-1
      FJ=J/XN
      DO 150 I=1,NR
      IF (WT(I) .EQ. 0.) GO TO 150
      KI=KEMP(I)
      WT(I)=WT(I)*FJ/TXN(KI)
  150 CONTINUE
      RETURN
      END

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


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,101
            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,error)
    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), error
    real(kind=8) :: start_time, end_time

!       call date_and_time(values=values_start)

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

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

!$OMP PARALLEL DO PRIVATE(I,J,K,sum,aamax,imax,d,dum) SHARED (A,indx)
        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
!$OMP END PARALLEL DO

!       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

! ---------------------------------------------------------
! This subroutine solves a linear system of equations using the SLAP routines
! It uses an Incomplete LU BiConjugate Gradient Sparse Ax=b solver
! ---------------------------------------------------------
!   subroutine slapsolver(a,b)
!   real(kind=8), INTENT(INOUT) :: a(:,:), b(:)
!   integer :: i, j, k, j2
!   real(kind=8), allocatable :: a_sparse(:), x(:), rwork(:)
!   integer, allocatable :: ia(:), ja(:), iwork(:)
!   integer :: n, nelt, isym, itol, itmax, iter, ierr, iunit, lenw, leniw, nsave
!   real(kind=8) :: tol, err
!   integer, allocatable :: indx(:)
!
!       n = size(b)
!       allocate(x(n))
!       k = 0
!       do i = 1, n
!           do j = 1, n
!               if (a(j,i) /= 0.d0) k = k + 1
!           enddo
!       enddo
!       nelt = k
!       if (verbose_mode == 1) then
!           print *, 'Number of non-zero elements : ', nelt
!       endif
!       allocate(a_sparse(nelt))
!       allocate(ia(nelt))
!       allocate(ja(nelt))
!
!       k = 1
!       do i = 1, n
!           do j = 1, n
!               if (a(j,i) /= 0.d0) then
!                   a_sparse(k) = a(j,i)
!                   ia(k) = j
!                   ja(k) = i
!                   k = k + 1
!               endif
!           enddo
!       enddo
!
!       x = 0.d0
!
!       isym = 0
!       itol = 1
!       tol = 1.d-10
!       itmax = 10000
!       iunit = 0
!       lenw = n*n
!       leniw = n*n
!       nsave = 10
!       allocate(iwork(leniw))
!       allocate(rwork(lenw))
!
!       if (verbose_mode == 1) then
!           print *, 'Solving the system...'
!       endif
!       call ds2y( n, nelt, ia, ja, a_sparse, isym )
!       call dslubc(n, b, x, nelt, ia, ja, a_sparse, isym, itol, tol, &
!           itmax, iter, err, ierr, iunit, rwork, lenw, iwork, leniw )
!       if (verbose_mode == 1) then
!           print *, 'Number of iterations : ', iter, 'of ', itmax
!           print *, 'Error : ', err, ierr
!       endif
!       if (iter >= itmax+1 .or. ierr /= 0) then
!           if (verbose_mode == 1) then
!               print *, 'Couldn''t reach convergence. Trying other initialization...'
!           endif
!           write(17,*) 'Couldn''t reach convergence. Trying other initialization...'
!           x = 1.d0
!           isym = 0
!           itol = 0
!           tol = 1.d-10
!           itmax = 10000
!           iunit = 0
!           lenw = n*n
!           leniw = n*n
!           nsave = 10
!           call dslucs(n, b, x, nelt, ia, ja, a_sparse, isym, itol, tol, &
!           itmax, iter, err, ierr, iunit, rwork, lenw, iwork, leniw )
!           if (verbose_mode == 1) then
!               print *, 'Number of iterations : ', iter, 'of ', itmax
!               print *, 'Error : ', err, ierr
!           endif
!       endif
!
!       if (iter >= itmax+1 .or. ierr /= 0) then
!           if (verbose_mode == 1) then
!               print *, 'Couldn''t reach convergence. Trying GMRES...'
!           endif
!           write(17,*) 'Couldn''t reach convergence. Trying GMRES...'
!           x = 0.d0
!           isym = 0
!           itol = 0
!           tol = 1.d-10
!           itmax = 10000
!           iunit = 0
!           lenw = n*n
!           leniw = n*n
!           nsave = 10
!           call dsdgmr(n, b, x, nelt, ia, ja, a_sparse, isym, nsave, itol, tol, &
!           itmax, iter, err, ierr, iunit, rwork, lenw, iwork, leniw )
!           if (verbose_mode == 1) then
!               print *, 'Number of iterations : ', iter, 'of ', itmax
!               print *, 'Error : ', err, ierr
!           endif
!       endif
!
!       if (iter >= itmax+1 .or. ierr /= 0) then
!           if (verbose_mode == 1) then
!               print *, 'Couldn''t reach convergence. Trying other initialization...'
!           endif
!           write(17,*) 'Couldn''t reach convergence. Trying other initialization...'
!           x = 1.d0
!           isym = 0
!           itol = 0
!           tol = 1.d-10
!           itmax = 10000
!           iunit = 0
!           lenw = n*n
!           leniw = n*n
!           nsave = 10
!           call dsdgmr(n, b, x, nelt, ia, ja, a_sparse, isym, nsave, itol, tol, &
!           itmax, iter, err, ierr, iunit, rwork, lenw, iwork, leniw )
!           if (verbose_mode == 1) then
!               print *, 'Number of iterations : ', iter, 'of ', itmax
!               print *, 'Error : ', err, ierr
!           endif
!       endif
!
!       b = x
!
!       deallocate(x)
!       deallocate(a_sparse)
!       deallocate(ia)
!       deallocate(ja)
!       deallocate(iwork)
!       deallocate(rwork)
!
!       if (iter >= itmax+1 .or. ierr /= 0) then
!           if (verbose_mode == 1) then
!               print *, 'Couldn''t reach convergence. Trying LU decomposition...'
!           endif
!           write(17,*) 'Couldn''t reach convergence. Trying LU decomposition...'
!           allocate(indx(n))
!           call ludcmp(a,indx,tol)
!           call lubksb(a,indx,b)
!           deallocate(indx)
!       endif
!
!   end subroutine slapsolver

! ---------------------------------------------------------
! 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

!--------------------------------------------------------------
! Dawson's integral
!--------------------------------------------------------------
    function dawson(dv)
    real(kind=8) :: dv(:)
    real(kind=8) :: dawson(size(dv)), c(6), x
    real(kind=8), parameter :: H=0.4, A1=2./3., A2=0.4, A3=2./7.
    integer :: i, n, j, n0
    real(kind=8) :: xx, xp, e1, e2, d1, d2, sum, x2

        n = size(dv)
        do i = 1, 6
            c(i)=exp(-((2.*float(i)-1.)*H)**2)
        enddo

        do i = 1, n
            if (abs(dv(i)) < 0.2) then
                x = dv(i)
                x2 = dv(i)**2
                dawson(i) = x*(1.-A1*x2*(1.-A2*x2*(1.-A3*x2)))
            else
                x = dv(i)
                xx = abs(dv(i))
                n0 = 2*nint(0.5*xx/H)
                xp = xx-float(n0)*H
                e1 = exp(2.*xp*H)
                e2 = e1**2
                d1 = float(n0+1)
                d2 = d1-2.
                sum = 0.
                do j = 1,6
                  sum = sum + c(j)*(e1/d1+1./(d2*e1))
                  d1 = d1+2.
                  d2 = d2-2.
                  e1 = e2*e1
                enddo
                dawson(i) = 0.5641895835*sign(exp(-xp**2),x)*sum
            endif
        enddo

    end function dawson

!--------------------------------------------------------------
! 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)

! For small a, evaluate the power expansion up to second order (Landi degl'Innocenti & Landolfi 2004, p. 168)
        if (da < 1.d-3) then
            profile = cmplx(exp(-dv**2) + 2.d0*da/sqrt(PI)*(2*dv*dawson(dv)-1.d0), 2.d0*dawson(dv) / sqrt(PI) - 2.d0*da*dv*exp(-dv**2))
            return
        endif
        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

        function erf(x)
        real(kind=8) :: erf, x, dumerfc, t, z

            z = abs(x)
            t = 1.0 / ( 1.0 + 0.5 * z )

            dumerfc =       t * exp(-z * z - 1.26551223 + t *       &
             ( 1.00002368 + t * ( 0.37409196 + t *      &
             ( 0.09678418 + t * (-0.18628806 + t *      &
             ( 0.27886807 + t * (-1.13520398 + t *      &
             ( 1.48851587 + t * (-0.82215223 + t * 0.17087277 )))))))))

            if ( x.lt.0.0 ) dumerfc = 2.0 - dumerfc

            erf = 1.0 - dumerfc

        end function erf

! ---------------------------------------------------------
! Solve a one-dimensional equality using the secant method
! ---------------------------------------------------------
    function confidenceLevel(x, nu, P)
    real(kind=8) :: confidenceLevel, x, nu, P
        confidenceLevel = gammq(0.5*nu,0.5*x) + P - 1.d0
        return
    end function confidenceLevel

! ---------------------------------------------------------
! Solve a one-dimensional equality using the secant method
! ---------------------------------------------------------
    function secantConfidenceLevel(nu,P)
    real(kind=8) :: func, x1, x2, xacc, fl, f, secantConfidenceLevel, xl, swap, dx, nu,P
    integer :: j
    integer, parameter :: maxit=30
        xacc = 1e-3
        x1 = nu-0.1
        x2 = nu+14.d0
        fl=confidenceLevel(x1,nu,P)
        f=confidenceLevel(x2,nu,P)
      if(abs(fl).lt.abs(f))then
        secantConfidenceLevel=x1
        xl=x2
        swap=fl
        fl=f
        f=swap
      else
        xl=x1
        secantConfidenceLevel=x2
      endif
      do j=1,maxit
        dx=(xl-secantConfidenceLevel)*f/(f-fl)
        xl=secantConfidenceLevel
        fl=f
        secantConfidenceLevel=secantConfidenceLevel+dx
        f=confidenceLevel(secantConfidenceLevel,nu,P)
        if(abs(dx).lt.xacc.or.f.eq.0.)return
      enddo
      print *, 'rtsec exceed maximum iterations'
      end function secantConfidenceLevel

! ---------------------------------------------------------
! Return the log of the gamma function
! ---------------------------------------------------------
      function gammln(xx)
      real(kind=8) :: cof(6),stp,half,one,fpf,x,tmp,ser,gammln, dx, xx
      data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,-1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
      data half,one,fpf/0.5d0,1.0d0,5.5d0/
      integer :: j
      x=xx-one
      tmp=x+fpf
      tmp=(x+half)*log(tmp)-tmp
      ser=one
      do j=1,6
        x=x+one
        ser=ser+cof(j)/x
      enddo
      gammln=tmp+log(stp*ser)
      return
      end function gammln

! ---------------------------------------------------------
!
! ---------------------------------------------------------
     subroutine gcf(gammcf,a,x,gln)
     real(kind=8) :: gammcf, a, x, gln, gold, a0, a1, b0, b1, fac, an, ana, anf, g
     real(kind=8), parameter :: eps=3e-7
     integer, parameter :: itmax=200
     integer :: n

      gln=gammln(a)
      gold=0.
      a0=1.
      a1=x
      b0=0.
      b1=1.
      fac=1.
      do n=1,itmax
        an=float(n)
        ana=an-a
        a0=(a1+a0*ana)*fac
        b0=(b1+b0*ana)*fac
        anf=an*fac
        a1=x*a0+anf*a1
        b1=x*b0+anf*b1
        if(a1.ne.0.)then
          fac=1./a1
          g=b1*fac
          if(abs((g-gold)/g).lt.eps) then
            gammcf=exp(-x+a*log(x)-gln)*g
            return
          endif
          gold=g
        endif
      enddo
        print *, 'a too large, itmax too small'
      end subroutine gcf

      subroutine gser(gamser,a,x,gln)
      real(kind=8) :: gamser, a, x, gln, ap, sum, del
      real(kind=8), parameter :: eps=3e-7
      integer, parameter :: itmax=200
      integer :: n

      gln=gammln(a)
      if(x.le.0.)then
        gamser=0.
        return
      endif
      ap=a
      sum=1./a
      del=sum
      do n=1,itmax
        ap=ap+1.
        del=del*x/ap
        sum=sum+del
        if(abs(del).lt.abs(sum)*eps) then
            gamser=sum*exp(-x+a*log(x)-gln)
            return
        endif
      enddo
      print *, 'a too large, itmax too small'
      end subroutine gser

    function gammq(a,x)
    real(kind=8) :: gammq, a, x, gln, gamser, gammcf
      if(x.lt.0..or.a.le.0.) then
        print *, x, a
        stop
    endif
      if(x.lt.a+1.)then
        call gser(gamser,a,x,gln)
        gammq=1.-gamser
      else
        call gcf(gammcf,a,x,gln)
        gammq=gammcf
      endif
      return
      end function gammq

end module maths

max_stars_repo_name: aasensio/hazel

max_issues_repo_name: aasensio/hazel

max_forks_repo_name: aasensio/hazel


text:
C    @(#)change.f	20.7 3/29/99
C****************************************************************
C
C   File: change.f
C
C   Purpose: Process /CHANGE commands
C
C   Author: BPA staff       Date: 18 May 1992
C   Called by: p_gtdata.f, prodat.f
C
C****************************************************************
C
      subroutine change ()

      include 'ipfinc/parametr.inc'

      include 'ipfinc/basval.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/changr.inc'
      include 'ipfinc/dtaiop.inc'
      include 'ipfinc/filnam.inc'
      include 'ipfinc/jobctl.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/oldchg.inc'
      include 'ipfinc/ordsta.inc'
      include 'ipfinc/prt.inc'
C
C       Update old base values
C
        call xdate(dte)
        basval(5) = dte
        call n_time( ihr, imin, isec )
        write (basval(6), 90) ihr, imin, isec
   90   format (i2, ':', i2, ':', i2)
        basval(8) = prgvsn
        basval(9) = usrnam

C
C       Set up headers
C
        call forbtm ()
        write (outbuf,100)
  100   format('* * * Network changes * * *')
        call rpnlod ()
        call fortop ()
        call space(2)
        call chread (numfil)

        if (numchg .gt. 0)   then
C
C          Convert to external order if in internal order.
C
           if (ordvlt .eq. 2) then
              ordvlt = 1
              call mvnew1d(e,opt2inp,ntot)
              call mvnew1d(f,opt2inp,ntot)
           endif
           if (ordcap .eq. 2) then
              ordcap = 1
              call mvnew2d(capcor,opt2inp,ntot)
           endif
           call chchek ()
           call chpct ()
           call chdel ()
           call chprea ()
           call chadd ()
           call chmods ()
           if (numchg .gt. 0) call sortbus ()
           call chmisc ()
           call chg_error ()
C
C          Format of oldchg()
C
C             (1:120) - change record
C
C           (121:121) - origination code
C                      ' ' indicates original record
C                      'P' inidcates pseudo-change record
C                      'T' indicates transpose record
C                      'U' indicates transpose psuedo-change record
c
C           (122:125) - chgcrd() record numnber, converted to an
c                       absolute index
C
C           (126:126) - process code
C                      ' ' indicates not processed
C                      'P' indicates processed without errors
C                      'E' indicates processed with errors
C
C           (127:128) - source of record code
C                      'n' indicates from imported change file no. "n"
C                      ' ' indicates interactive change
C
           do ic = 1, numchg
              call storchgs (chgcrd(ic), numfil)
           enddo
        else
           outbuf = '0 No change records encountered'
           call prtout(1)
        endif

        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