text:
!! This module contains tools to convert fortran object
module mod_convert
  !! This module contains tools to convert fortran object
  implicit none

  private

  public :: str

contains

  character(len=20) function str(k)
    !! Convert an integer to string
    integer, intent(in) :: k
    write (str, *) k
    str = adjustl(str)
  end function str

end module mod_convert

max_stars_repo_name: antoinemarchal/ROHSA

max_issues_repo_name: antoinemarchal/ROHSA

max_forks_repo_name: antoinemarchal/ROHSA


text:
program test
   use numz
   use ran_mod
   real(b8)x
   integer iseed
   iseed=-1234
   x=ran1(iseed)
   do i=1,10
      write(*,*)i,ran1()
   enddo
end program

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
Viet is known to mooch food from the FAIS room.

Viet likes stuffed animals. These animals are named Elephant Elephant and bee bees

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!! main subroutines for mt measurement
!

! =================================================================================================
! subroutine mt_measure()
! Boxcar/Cosine/Multitaper estimates of the phase function between dat1 and dat2 (dat1-dat2)
!
!  Input:
!        dat1_dt(:), dat2_dt(:), t0, dt, npts -- original dat1 and dat2 array
!        tstart, tend -- start and end of the measurement window (can be from Alessia's code)
!  Output:
!        istart -- starting index of the windowed portion of  original trace
!        dat1_dtw(:), dat2_dtw(:), nlen -- windowed and shifted dat1, windowed dat2
!        tshift, dlnA, cc_max -- time shift and amplitude cross-correlation measurements
!        i_right -- the maximum reliable frequency estimate index
!        dtau_w(:), dlnA_w(:) -- estimates of travel-time and amplitude anomaly
!        err_dt(:), err_dlnA(:) -- error bar of the travel-time and amplitude estimates (MT only)
!
! =================================================================================================

subroutine mt_measure(dat1,dat2,npts,deltat,nlen,tshift_cc,dlnA_cc, i_fstart,i_fend,&
        wvec,&
        trans_func,dtau_w,dlnA_w,err_dtau,err_dlnA)
    use constants
    implicit none

    integer, intent(in) :: nlen, npts
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: dat1,dat2
    real(kind=CUSTOM_REAL), intent(in) ::  deltat,tshift_cc,dlnA_cc
    integer, intent(in) :: i_fstart,i_fend
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: wvec
    complex (CUSTOM_REAL), dimension(NPT), intent(out) :: trans_func
    real(kind=CUSTOM_REAL), dimension(*), intent(out) :: dtau_w,dlnA_w,err_dtau,err_dlnA

    !! fft -- double size
    ! FFT parameters
    real(kind=SIZE_DOUBLE) :: ampmax,wtr_mtm_use,wtr_use
    integer :: i,ictaper,iom

    !! multitaper
    real(kind=SIZE_DOUBLE), dimension(NPT) :: ey1,ey2
    real(kind=SIZE_DOUBLE), dimension(:,:),allocatable :: tas

    !! fft to estimate transfer function
    ! tapered data in time domain
    real(kind=CUSTOM_REAL), dimension(NPT) :: dat1_dtw_h,dat2_dtw_h
    ! fft domain
    complex (SIZE_DOUBLE), dimension(NPT) :: dat2_dtwo, dat1_dtwo, &
        dat2_dtw_ho, dat1_dtw_ho,  &
        top_mtm,bot_mtm
    !! transfer function to phase/amplitude measurement
    real(kind=CUSTOM_REAL), dimension(NPT) :: phi_mtm,dtau_mtm,abs_mtm,dlnA_mtm
    real(kind=CUSTOM_REAL), dimension(:,:),allocatable :: phi_mul,dtau_mul,abs_mul,dlnA_mul

    !! jackknife error
    real(kind=CUSTOM_REAL) :: edt_ave,edt_iom,edA_ave,edA_iom

    character (len=MAX_FILENAME_LEN) :: filename

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

    if ( nlen > npts .or. nlen <1) then
        print*, nlen, npts
        stop 'Check nlen'
    endif

    !-------------------------------------------------------------------------------
    ! multitaper estimation of transfer function
    !-------------------------------------------------------------------------------
    !print*,'nlen,NW,ntaper -- ',nlen, NW, ntaper
    ! calculate the tapers
    allocate(tas(NPT,ntaper))
    call staper(nlen, dble(NW), ntaper, tas, NPT, ey1, ey2)

    if( DISPLAY_DETAILS) then
        do ictaper=1,ntaper
        write(filename,'(i)') ictaper
        open(1,file=trim(output_dir)//'/'//'taper_t_'//trim(adjustl(filename)),status='unknown')
        open(2,file=trim(output_dir)//'/'//'taper_w_'//trim(adjustl(filename)),status='unknown')
        do i = 1,nlen
        write(1,'(f15.2,E15.5)') (i-1)*deltat,tas(i,ictaper)
        ! write(1,'(f15.2,<ntaper>e15.5)') (i-1)*deltat,tas(i,1:ntaper)
        enddo
        close(1)
        dat1_dtw_ho(:) = cmplx(0.0_SIZE_DOUBLE,0.0_SIZE_DOUBLE)
        dat1_dtw_h(1:nlen) = real(tas(1:nlen,ictaper),CUSTOM_REAL)
        dat1_dtw_ho(1:nlen) = dcmplx(dat1_dtw_h(1:nlen),0.0_SIZE_DOUBLE)
        call fft(LNPT,dat1_dtw_ho,dble(FORWARD_FFT),dble(deltat))
        do  i =  1,50
        write(2,'(2f15.2,E15.5)') wvec(i),abs(dat1_dtw_ho(i))
        enddo
        close(1)
        close(2)
        enddo
    endif

    !! estimate transfer function from mtm
    ! initialize transfer function terms
    top_mtm(:)   = cmplx(0.0_SIZE_DOUBLE,0.0_SIZE_DOUBLE)
    bot_mtm(:)   = cmplx(0.0_SIZE_DOUBLE,0.0_SIZE_DOUBLE)
    trans_func(:) = cmplx(0.0_SIZE_DOUBLE,0.0_SIZE_DOUBLE)

    do ictaper = 1, ntaper

    dat2_dtw_ho(:) = cmplx(0.0_SIZE_DOUBLE,0.0_SIZE_DOUBLE) ! note: this has to be initialized inside the loop
    dat1_dtw_ho(:) = cmplx(0.0_SIZE_DOUBLE,0.0_SIZE_DOUBLE)

    ! apply time-domain taper
    do i = 1, nlen
    dat2_dtw_h(i) = dat2(i) * real(tas(i,ictaper),CUSTOM_REAL)     ! single-tapered, windowed data2
    dat1_dtw_h(i) = dat1(i) * real(tas(i,ictaper),CUSTOM_REAL)  ! single-tapered, windowed, shifted data1
    enddo

    dat2_dtw_ho(1:nlen) = dcmplx(dat2_dtw_h(1:nlen),0.0_SIZE_DOUBLE)
    dat1_dtw_ho(1:nlen) = dcmplx(dat1_dtw_h(1:nlen),0.0_SIZE_DOUBLE)

    ! apply FFT to get complex spectra
    call fft(LNPT,dat2_dtw_ho,dble(FORWARD_FFT),dble(deltat))
    call fft(LNPT,dat1_dtw_ho,dble(FORWARD_FFT),dble(deltat))

    if(DISPLAY_DETAILS .and. ictaper==1) then
        open(1,file=trim(output_dir)//'/fft_in',status='unknown')
        open(2,file=trim(output_dir)//'/fft_out',status='unknown')
        do i = 1, nlen
        write(1,'(3E15.5)') i,dat1_dtw_h(i),dat2_dtw_h(i)
        enddo
        do  i =  i_fstart,i_fend
        write(2,'(3E15.5)') wvec(i)/TWOPI, abs(dat1_dtw_ho(i)),abs(dat2_dtw_ho(i))
        enddo
        close(1)
        close(2)
    endif

    ! calculate top and bottom of transfer function
    do i = i_fstart,i_fend
    top_mtm(i) = top_mtm(i) + dat1_dtw_ho(i) * conjg(dat2_dtw_ho(i))
    bot_mtm(i) = bot_mtm(i) + dat2_dtw_ho(i) * conjg(dat2_dtw_ho(i))
    enddo

    enddo  ! ictapers

    ! water level for transfer function
    ampmax = maxval(abs(bot_mtm))

    ! LQY -- is this too small ???
    wtr_mtm_use = ampmax * wtr_mtm**2
    wtr_mtm_use = ampmax * 0.01

    ! calculate MT transfer function using water level
    do i =  i_fstart,i_fend
    !      if(abs(bot_mtm(i)) > abs(wtr_mtm_use)) trans_func(i) = cmplx(top_mtm(i) / bot_mtm(i))
    !      if(abs(bot_mtm(i)) < abs(wtr_mtm_use)) trans_func(i) = cmplx(top_mtm(i) / (bot_mtm(i)+wtr_mtm_use))
    trans_func(i) = cmplx(top_mtm(i) / (bot_mtm(i)+wtr_mtm_use))
    enddo

    if(DISPLAY_DETAILS) then
        ! print*
        ! print*,'water level for multi taper transfer function is : ',wtr_use
        ! fft of untapered inputs
        dat1_dtwo = cmplx(0.0_SIZE_DOUBLE,0.0_SIZE_DOUBLE)
        dat2_dtwo = cmplx(0.0_SIZE_DOUBLE,0.0_SIZE_DOUBLE)
        dat1_dtwo(1:nlen) = dcmplx(dat1(1:nlen),0.0)
        dat2_dtwo(1:nlen) = dcmplx(dat2(1:nlen),0.0)
        call fft(LNPT,dat1_dtwo,dble(FORWARD_FFT),dble(deltat))
        call fft(LNPT,dat2_dtwo,dble(FORWARD_FFT),dble(deltat))
        open(1,file=trim(output_dir)//'/top_bot_mtm',status='unknown')
        open(2,file=trim(output_dir)//'/trans_func',status='unknown')
        open(3,file=trim(output_dir)//'/trans_error',status='unknown')
        do  i =  i_fstart,i_fend
        write(1,'(3E15.5)') wvec(i)/TWOPI, abs(top_mtm(i)),abs(bot_mtm(i))
        write(2,'(2E15.5)') wvec(i)/TWOPI, abs(trans_func(i))
        write(3,'(4E15.5)') wvec(i)/TWOPI, abs(dat1_dtwo(i)),abs(dat2_dtwo(i)),abs(dcmplx(trans_func(i))*dat2_dtwo(i))
        enddo
        close(1)
        close(2)
        close(3)
    endif

    !! phase and amplitude measurement derived from transfer functions
    phi_mtm(:) = 0.d0
    abs_mtm(:) = 0.d0
    dtau_mtm(:) = 0.d0
    dlnA_mtm(:) = 0.d0
    call write_trans(trans_func, wvec, i_fstart,i_fend,tshift_cc,dlnA_cc,phi_mtm, abs_mtm,dtau_mtm, dlnA_mtm)

    ! pass results to main routine
    dtau_w(i_fstart:i_fend) = dtau_mtm(i_fstart:i_fend)
    dlnA_w(i_fstart:i_fend) = dlnA_mtm(i_fstart:i_fend)

    if(DISPLAY_DETAILS) then
        open(1,file=trim(output_dir)//'/dtau_dlnA',status='unknown')
        do  i =  i_fstart,i_fend
        write(1,'(3E15.5)') wvec(i)/TWOPI,dtau_w(i),dlnA_w(i)
        enddo
        close(1)
    endif

    !-------------------------------------------------------------------------------
    ! multitaper error estimation
    !-------------------------------------------------------------------------------

    if (ntaper > 1 .and. USE_ERROR_MT) then

        ! allocate Jacknife MT estimates
        allocate(phi_mul(NPT,ntaper))
        allocate(dtau_mul(NPT,ntaper))
        allocate(abs_mul(NPT,ntaper))
        allocate(dlnA_mul(NPT,ntaper))

        do iom = 1, ntaper

        top_mtm(:) = cmplx(0.,0.)
        bot_mtm(:) = cmplx(0.,0.)

        do ictaper = 1, ntaper
        if(ictaper.eq.iom) cycle

        ! apply ictaper-th taper
        dat2_dtw_h(1:nlen) = dat2(1:nlen) * real(tas(1:nlen,ictaper),CUSTOM_REAL)
        dat1_dtw_h(1:nlen) = dat1(1:nlen) * real(tas(1:nlen,ictaper),CUSTOM_REAL)

        ! complex tapered series
        dat2_dtw_ho(:) = cmplx(0.,0.)
        dat1_dtw_ho(:) = cmplx(0.,0.)
        dat2_dtw_ho(1:nlen) = dcmplx(dat2_dtw_h(1:nlen),0.)
        dat1_dtw_ho(1:nlen) = dcmplx(dat1_dtw_h(1:nlen),0.)

        ! apply f.t. to get complex spectra
        call fft(LNPT,dat2_dtw_ho,dble(FORWARD_FFT),dble(deltat))
        call fft(LNPT,dat1_dtw_ho,dble(FORWARD_FFT),dble(deltat))

        ! calculate top and bottom of Jacknife transfer function
        do i = i_fstart,i_fend
        top_mtm(i) = top_mtm(i) + dat1_dtw_ho(i) * conjg(dat2_dtw_ho(i))
        bot_mtm(i) = bot_mtm(i) + dat2_dtw_ho(i) * conjg(dat2_dtw_ho(i))
        enddo
        enddo ! ictaper

        ! water level
        ampmax = maxval(abs(bot_mtm))
        wtr_use = cmplx(ampmax * wtr_mtm ** 2, 0.)
        wtr_mtm_use = ampmax * 0.01

        !  calculate transfer function using water level
        do i = i_fstart,i_fend
        !if(abs(bot_mtm(i)).gt.abs(wtr_use)) trans_func(i) = top_mtm(i) / bot_mtm(i)
        !if(abs(bot_mtm(i)).le.abs(wtr_use)) trans_func(i) = top_mtm(i) /(bot_mtm(i)+wtr_use)
        trans_func(i) = cmplx(top_mtm(i) /(bot_mtm(i)+wtr_use))
        enddo

        call write_trans(trans_func, wvec,i_fstart,i_fend,tshift_cc,dlnA_cc,&
            phi_mul(:,iom),abs_mul(:,iom),dtau_mul(:,iom), dlnA_mul(:,iom))

        enddo ! iom

        !----------------------
        err_dtau(1:NPT)   = 0.
        err_dlnA(1:NPT)   = 0.

        do i = i_fstart,i_fend

        edt_ave   = 0.
        edA_ave   = 0.

        do iom = 1, ntaper
        edt_iom = ntaper*dtau_mtm(i) - (ntaper-1)*dtau_mul(i,iom)
        edt_ave = edt_ave + edt_iom

        edA_iom = ntaper*dlnA_mtm(i) - (ntaper-1)*dlnA_mul(i,iom)
        edA_ave = edA_ave + edA_iom

        enddo

        edt_ave   = edt_ave   / (ntaper)
        edA_ave   = edA_ave   / (ntaper)

        do iom = 1, ntaper
        err_dlnA(i)  = err_dlnA(i) + ( dlnA_mul(i,iom) - edA_ave)**2
        err_dtau(i)   = err_dtau(i)  + (dtau_mul(i,iom) - edt_ave)**2
        enddo

        err_dlnA(i)  =  sqrt( err_dlnA(i) / (ntaper * (ntaper-1) ) )
        err_dtau(i)   =  sqrt( err_dtau(i) / (ntaper * (ntaper-1) ) )
        ! set the error bar for the first point corresponding to
        ! static offset to be large, which makes no contribution to
        ! the adjoint source
        !  if (i == 1) err_dt(i) = LARGE_VAL

        enddo ! i_fstart

        deallocate(phi_mul)
        deallocate(dtau_mul)
        deallocate(abs_mul)
        deallocate(dlnA_mul)

        if(DISPLAY_DETAILS) then
            open(1,file=trim(output_dir)//'/err_dtau_dlnA',status='unknown')
            do  i =  i_fstart,i_fend
            write(1,'(3E15.5)') wvec(i)/TWOPI,err_dtau(i),err_dlnA(i)
            enddo
            close(1)
        endif

    endif  ! if use_error_MT

    !     ------------------------------------------------------------------
    !     End error calculation loop
    !     ------------------------------------------------------------------

    deallocate(tas)

end subroutine mt_measure
! ------------------------------------------------------------------------------------
subroutine mtm_adj(syn,npts,deltat,nlen,df,i_fstart,i_fend,dtau_w,dlnA_w,&
        err_dt_cc,err_dlnA_cc, &
        err_dtau_mt,err_dlnA_mt, &
        fp,fq)
    use constants
    implicit none

    integer, intent(in) :: npts,nlen
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: syn
    real(kind=CUSTOM_REAL), intent(in) ::  deltat,df,err_dt_cc,err_dlnA_cc
    integer, intent(in) :: i_fstart,i_fend
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: dtau_w,dlnA_w,err_dtau_mt,err_dlnA_mt
    real(kind=CUSTOM_REAL), dimension(*), intent(out) :: fp, fq

    ! frequency-domain taper
    integer :: nfrange, window_type =3
    real(kind=CUSTOM_REAL), dimension(NPT) :: w_taper, wp_taper, wq_taper
    real(kind=CUSTOM_REAL) :: ffac, dtau_wtr, dlnA_wtr, err_t, err_A

    !! fft -- double size
    ! FFT parameters
    real(kind=SIZE_DOUBLE) :: ampmax,wtr_mtm_use,wtr_use
    integer :: i,ictaper,iom

    !! multitaper
    real(kind=SIZE_DOUBLE), dimension(NPT) :: ey1,ey2
    real(kind=SIZE_DOUBLE), dimension(:,:),allocatable :: tas

    !! fft to evaluate adj
    ! tapered data in time domain
    real(kind=CUSTOM_REAL), dimension(NPT) :: syn_dtw_h, syn_vtw_h
    ! fft and ifft (forward and backward)
    complex (SIZE_DOUBLE), dimension(:,:),allocatable :: syn_dtw_ho_all,syn_vtw_ho_all
    complex (SIZE_DOUBLE), dimension(NPT) :: p_bot_mtm, q_bot_mtm, &
        pwc_adj,qwc_adj
    real(kind=SIZE_DOUBLE), dimension(NPT) ::  dtau_pj_t,dlnA_qj_t

    ! frequency-domain tapers
    ! THIS CHOICE WILL HAVE AN EFFECT ON THE ADJOINT SOURCES
    nfrange = i_fend -i_fstart
    w_taper(:) = 0.
    do i = i_fstart, i_fend
    ! if(window_type ==1) w_taper(i) = 1. ! boxcar window
    ! if(window_type ==2) w_taper(i) = 1. - (2.0/nw)**2 * ((i-1) - nw/2.0)**2       !welch window
    if(window_type ==3) w_taper(i) = &
        1. - cos(PI*(i-i_fstart)/(nfrange))**ipwr_w !cosine window
    enddo

    ! compute normalization factor for w_taper
    ! note: 2 is needed for the integration from -inf to inf
    ffac = 2.0 * df * sum(w_taper(i_fstart:i_fend) )
    if (DISPLAY_DETAILS) then
        print*
        print*, ' normalize frequency-domain taper'
        print*, ' df = ', df
        print *, 'Taper normalization factor, ffac = ', ffac
    endif
    ! no error estimate
    ! only adds normalization factor
    wp_taper = w_taper / ffac
    wq_taper = w_taper / ffac

    ! add error estimated
    if (USE_ERROR_CC) then
        !          ! CC error estimate
        do i = i_fstart,i_fend
        wp_taper(i) = wp_taper(i) / (err_dt_cc ** 2)
        wq_taper(i) = wq_taper(i) / (err_dlnA_cc ** 2)
        enddo
    endif
    if (USE_ERROR_MT) then
        ! MT jack-knife error estimate
        dtau_wtr = WTR * sum(abs(dtau_w(i_fstart:i_fend)))/ nfrange
        dlnA_wtr = WTR * sum(abs(dlnA_w(i_fstart:i_fend)))/ nfrange
        print*
        print*,'to add mt error :'
        print*,'water level for err_dtau is ',dtau_wtr
        print*,'water level for err_dlnA is ',dlnA_wtr

        do i = i_fstart,i_fend
        err_t = err_dtau_mt(i)
        err_A = err_dlnA_mt(i)
        if (err_t < dtau_wtr)  err_t = err_t + dtau_wtr
        if (err_A < dlnA_wtr)  err_A = err_A + dlnA_wtr
        wp_taper(i) = wp_taper(i) / (err_t ** 2)
        wq_taper(i) = wq_taper(i) / (err_A ** 2)
        enddo
    endif

    if( DISPLAY_DETAILS) then
        open(1,file=trim(output_dir)//'/frequency_taper',status='unknown')
        do  i =  i_fstart,i_fend
        write(1,'(I5,3e15.5)') i,w_taper(i),wp_taper(i),wq_taper(i)
        enddo
        close(1)
    endif

    ! allocate MT variables
    allocate(syn_dtw_ho_all(NPT,ntaper))
    allocate(syn_vtw_ho_all(NPT,ntaper))

    ! calculate the tapers
    allocate(tas(NPT,ntaper))
    call staper(nlen, NW, ntaper, tas, NPT, ey1, ey2)

    p_bot_mtm = 0.
    q_bot_mtm = 0.

    do ictaper = 1,ntaper

    ! tapered synthetic displacement
    syn_dtw_h(1:nlen) = syn(1:nlen) * real(tas(1:nlen,ictaper),CUSTOM_REAL)

    ! compute velocity of tapered syn
    call compute_vel(syn_dtw_h,npts,deltat,nlen,syn_vtw_h)

    ! single-tapered complex synthetic displacement and velocity
    syn_dtw_ho_all(:,ictaper) = 0.
    syn_vtw_ho_all(:,ictaper) = 0.
    syn_dtw_ho_all(1:nlen,ictaper) = dcmplx(syn_dtw_h(1:nlen),0.)
    syn_vtw_ho_all(1:nlen,ictaper) = dcmplx(syn_vtw_h(1:nlen),0.)

    ! apply FFT get complex spectra
    call fft(LNPT,syn_dtw_ho_all(:,ictaper),dble(FORWARD_FFT),dble(deltat))
    call fft(LNPT,syn_vtw_ho_all(:,ictaper),dble(FORWARD_FFT),dble(deltat))

    p_bot_mtm(:) = p_bot_mtm(:) + syn_vtw_ho_all(:,ictaper) &
        * conjg(syn_vtw_ho_all(:,ictaper))
    q_bot_mtm(:) = q_bot_mtm(:) + syn_dtw_ho_all(:,ictaper) &
        * conjg(syn_dtw_ho_all(:,ictaper))

    enddo ! ictaper

    if( DISPLAY_DETAILS) then
        open(2,file=trim(output_dir)//'/adj_bot_pq',status='unknown')
        do  i = i_fstart,i_fend
        write(2,*) i,abs(p_bot_mtm(i)), abs(q_bot_mtm(i))
        enddo
        close(2)
    endif

    ampmax = maxval(abs(p_bot_mtm))
    !    wtr_use = ampmax * wtr_mtm**2
    wtr_use = ampmax * 0.01

    ! compute p_j, q_j, P_j, Q_j and adjoint source fp, fq
    fp(1:npts) = 0.
    fq(1:npts) = 0.
    do ictaper = 1,ntaper

    ! compute p_j(w) and q_j(w)
    pwc_adj(:) = cmplx(0.,0.)
    qwc_adj(:) = cmplx(0.,0.)

    do i = i_fstart,i_fend
    ! if(abs(p_bot_mtm(i)) > abs(wtr_use)) pwc_adj(i) = &
    !        syn_vtw_ho_all(i,ictaper) / p_bot_mtm(i)
    ! if(abs(p_bot_mtm(i)) > abs(wtr_use)) pwc_adj(i) = &
    !        syn_vtw_ho_all(i,ictaper) / (p_bot_mtm(i) + wtr_use)
    pwc_adj(i) =  syn_vtw_ho_all(i,ictaper) / p_bot_mtm(i)
    qwc_adj(i) = -syn_dtw_ho_all(i,ictaper) / q_bot_mtm(i)
    enddo

    ! compute P_j(w) and Q_j(w)
    ! NOTE: the MT measurement is incorporated here
    !             also note that wp_taper and wq_taper can contain
    !             uncertainty estimations
    ! adds misfit measurement dtau, dlnA
    pwc_adj(i_fstart: i_fend) = pwc_adj(i_fstart: i_fend) &
        * dcmplx(dtau_w(i_fstart: i_fend),0.) &
        * dcmplx(wp_taper(i_fstart: i_fend),0.)
    qwc_adj(i_fstart: i_fend) = qwc_adj(i_fstart: i_fend) &
        * dcmplx(dlnA_w(i_fstart: i_fend),0.) &
        * dcmplx(wq_taper(i_fstart: i_fend),0.)

    ! IFFT into the time domain
    call fftinv(LNPT,pwc_adj,dble(REVERSE_FFT),dble(deltat),dtau_pj_t)
    call fftinv(LNPT,qwc_adj,dble(REVERSE_FFT),dble(deltat),dlnA_qj_t)

    ! create adjoint source
    ! applies taper to time signal
    fp(1:nlen) = fp(1:nlen) + real(tas(1:nlen,ictaper) * dtau_pj_t(1:nlen),CUSTOM_REAL)
    fq(1:nlen) = fq(1:nlen) + real(tas(1:nlen,ictaper) * dlnA_qj_t(1:nlen),CUSTOM_REAL)

    enddo  ! ictaper
    deallocate(tas)
    deallocate(syn_dtw_ho_all)
    deallocate(syn_vtw_ho_all)

end subroutine mtm_adj
! ------------------------------------------------------------------------------------
subroutine mtm_DD_adj(s1,s2,npts,deltat,nlen,df,i_fstart,i_fend,ddtau_w,ddlnA_w,&
        err_dt_cc_obs,err_dt_cc_syn,err_dlnA_cc_obs,err_dlnA_cc_syn, &
        err_dtau_mt_obs,err_dtau_mt_syn,err_dlnA_mt_obs,err_dlnA_mt_syn, &
        fp1,fp2,fq1,fq2)
    use constants
    implicit none

    integer, intent(in) :: npts,nlen
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: s1,s2
    real(kind=CUSTOM_REAL), intent(in) :: deltat,df
    real(kind=CUSTOM_REAL), intent(in) :: err_dt_cc_obs,err_dt_cc_syn,err_dlnA_cc_obs,err_dlnA_cc_syn
    integer, intent(in) :: i_fstart,i_fend
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: ddtau_w,ddlnA_w
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: err_dtau_mt_obs,err_dtau_mt_syn,err_dlnA_mt_obs,err_dlnA_mt_syn
    real(kind=CUSTOM_REAL), dimension(*), intent(out) :: fp1,fp2,fq1,fq2

    ! frequency-domain taper
    integer :: nfrange, window_type =3
    real(kind=CUSTOM_REAL), dimension(NPT) :: w_taper, wp_taper, wq_taper
    real(kind=CUSTOM_REAL) :: ffac, ddtau_wtr, ddlnA_wtr, err_t, err_A

    !! fft -- double size
    ! FFT parameters
    real(kind=SIZE_DOUBLE) :: ampmax,wtr_mtm_use,wtr_use
    integer :: i,ictaper,iom

    !! multitaper
    real(kind=SIZE_DOUBLE), dimension(NPT) :: ey1,ey2
    real(kind=SIZE_DOUBLE), dimension(:,:),allocatable :: tas

    !! fft to evaluate adj
    ! tapered data in time domain
    real(kind=CUSTOM_REAL), dimension(NPT) :: s1_dtw_h, s1_vtw_h
    real(kind=CUSTOM_REAL), dimension(NPT) :: s2_dtw_h, s2_vtw_h
    ! fft and ifft (forward and backward)
    complex (SIZE_DOUBLE), dimension(:,:),allocatable :: s1_dtw_ho_all,s1_vtw_ho_all
    complex (SIZE_DOUBLE), dimension(:,:),allocatable :: s2_dtw_ho_all,s2_vtw_ho_all
    complex (SIZE_DOUBLE), dimension(NPT) :: Mtr1_mtm,Mtr2_mtm,Mtr3_mtm,Mtr4_mtm,Mtr5_mtm
    complex (SIZE_DOUBLE), dimension(NPT) :: p1_adj,pw1_adj,p2_adj,pw2_adj
    complex (SIZE_DOUBLE), dimension(NPT) :: q1_adj,qw1_adj,q2_adj,qw2_adj
    complex (SIZE_DOUBLE), dimension(NPT) :: fp1_adj,fp2_adj,fq1_adj,fq2_adj
    real(kind=SIZE_DOUBLE), dimension(NPT) :: fp1_adj_t,fp2_adj_t,fq1_adj_t,fq2_adj_t

    ! frequency-domain taper
    !    integer :: nfrange, window_type =3
    !    real(kind=CUSTOM_REAL), dimension(NPT) :: w_taper, wp_taper, wq_taper, ey1, ey2
    !    real(kind=CUSTOM_REAL) :: ffac, ddtau_wtr, ddlnA_wtr, err_t, err_A

    ! mtm adj
    !    real(kind=CUSTOM_REAL), dimension(NPT) :: s1_dtw_h, s1_vtw_h,s2_dtw_h, s2_vtw_h
    !    complex (CUSTOM_REAL), dimension(:,:),allocatable ::  s1_dtw_ho_all,s1_vtw_ho_all
    !    complex (CUSTOM_REAL), dimension(:,:),allocatable ::  s2_dtw_ho_all,s2_vtw_ho_all
    !    comple (CUSTOM_REAL), dimension(NPT) :: Mtr1_mtm,Mtr2_mtm,Mtr3_mtm,Mtr4_mtm,Mtr5_mtm
    !    complex (CUSTOM_REAL), dimension(NPT) :: p1_adj,pw1_adj,p2_adj,pw2_adj
    !    complex (CUSTOM_REAL), dimension(NPT) :: q1_adj,qw1_adj,q2_adj,qw2_adj
    !    complex (CUSTOM_REAL), dimension(NPT) :: fp1_adj,fp2_adj,fq1_adj,fq2_adj
    !    real(kind=CUSTOM_REAL), dimension(NPT) :: fp1_adj_t,fp2_adj_t,fq1_adj_t,fq2_adj_t
    !    real(kind=CUSTOM_REAL), dimension(:,:),allocatable :: tas

    ! frequency-domain tapers
    ! THIS CHOICE WILL HAVE AN EFFECT ON THE ADJOINT SOURCES
    nfrange = i_fend -i_fstart
    w_taper(:) = 0.
    do i = i_fstart, i_fend
    ! boxcar
    ! if(window_type ==1) w_taper(i) = 1. ! boxcar window
    ! welch
    ! if(window_type ==2) w_taper(i) = 1. - (2.0/nw)**2 * ((i-1) - nw/2.0)**2
    ! cosine
    if(window_type ==3) w_taper(i) = &
        1. - cos(PI*(i-i_fstart)/(nfrange))**ipwr_w
    enddo

    ! compute normalization factor for w_taper
    ! note: 2 is needed for the integration from -inf to inf
    ffac = 2.0 * df * sum(w_taper(i_fstart:i_fend) )
    if (DISPLAY_DETAILS) then
        print*
        print*, ' normalize frequency-domain taper'
        print*, ' df = ', df
        print *, 'Taper normalization factor, ffac = ', ffac
    endif
    ! no error estimate
    ! only adds normalization factor
    wp_taper = w_taper / ffac
    wq_taper = w_taper / ffac

    ! add error estimated
    if (USE_ERROR_CC) then
        !          ! CC error estimate
        do i = i_fstart,i_fend
        wp_taper(i) = wp_taper(i) / (err_dt_cc_obs * err_dt_cc_syn)
        wq_taper(i) = wq_taper(i) / (err_dlnA_cc_obs * err_dlnA_cc_syn)
        enddo
    endif
    if (USE_ERROR_MT) then
        ! MT jack-knife error estimate
        ddtau_wtr = (WTR * sum(abs(ddtau_w(i_fstart:i_fend)))/ nfrange)**2
        ddlnA_wtr = (WTR * sum(abs(ddlnA_w(i_fstart:i_fend)))/ nfrange)**2
        print*
        print*,'to add mt error :'
        print*,'water level for err_dtau is ',ddtau_wtr
        print*,'water level for err_dlnA is ',ddlnA_wtr

        do i = i_fstart,i_fend
        err_t = err_dtau_mt_obs(i) * err_dtau_mt_syn(i)
        err_A = err_dlnA_mt_obs(i) * err_dlnA_mt_syn(i)
        if (err_t < ddtau_wtr)  err_t = err_t + ddtau_wtr
        if (err_A < ddlnA_wtr)  err_A = err_A + ddlnA_wtr
        wp_taper(i) = wp_taper(i) / err_t
        wq_taper(i) = wq_taper(i) / err_A
        enddo
    endif

    if( DISPLAY_DETAILS) then
        open(1,file=trim(output_dir)//'/frequency_taper',status='unknown')
        do  i =  i_fstart,i_fend
        write(1,'(I5,3e15.5)') i,w_taper(i),wp_taper(i),wq_taper(i)
        enddo
        close(1)
    endif

    ! allocate MT variables
    allocate(s1_dtw_ho_all(NPT,ntaper))
    allocate(s1_vtw_ho_all(NPT,ntaper))
    allocate(s2_dtw_ho_all(NPT,ntaper))
    allocate(s2_vtw_ho_all(NPT,ntaper))

    ! calculate the tapers
    allocate(tas(NPT,ntaper))
    call staper(nlen, NW, ntaper, tas, NPT, ey1, ey2)

    !! constant terms for adj
    Mtr1_mtm = 0.
    Mtr2_mtm = 0.
    Mtr3_mtm = 0.
    Mtr4_mtm = 0.
    Mtr5_mtm = 0.

    do ictaper = 1,ntaper

    ! tapered synthetic displacement
    s1_dtw_h(1:nlen) = s1(1:nlen) * real(tas(1:nlen,ictaper),CUSTOM_REAL)
    s2_dtw_h(1:nlen) = s2(1:nlen) * real(tas(1:nlen,ictaper),CUSTOM_REAL)

    ! compute velocity of tapered syn
    call compute_vel(s1_dtw_h,npts,deltat,nlen,s1_vtw_h)
    call compute_vel(s2_dtw_h,npts,deltat,nlen,s2_vtw_h)

    ! single-tapered complex synthetic displacement and velocity
    s1_dtw_ho_all(:,ictaper) = 0.
    s1_vtw_ho_all(:,ictaper) = 0.
    s2_dtw_ho_all(:,ictaper) = 0.
    s2_vtw_ho_all(:,ictaper) = 0.
    s1_dtw_ho_all(1:nlen,ictaper) = dcmplx(s1_dtw_h(1:nlen),0.)
    s1_vtw_ho_all(1:nlen,ictaper) = dcmplx(s1_vtw_h(1:nlen),0.)
    s2_dtw_ho_all(1:nlen,ictaper) = dcmplx(s2_dtw_h(1:nlen),0.)
    s2_vtw_ho_all(1:nlen,ictaper) = dcmplx(s2_vtw_h(1:nlen),0.)

    ! apply FFT get complex spectra
    call fft(LNPT,s1_dtw_ho_all(:,ictaper),dble(FORWARD_FFT),dble(deltat))
    call fft(LNPT,s1_vtw_ho_all(:,ictaper),dble(FORWARD_FFT),dble(deltat))
    call fft(LNPT,s2_dtw_ho_all(:,ictaper),dble(FORWARD_FFT),dble(deltat))
    call fft(LNPT,s2_vtw_ho_all(:,ictaper),dble(FORWARD_FFT),dble(deltat))

    Mtr1_mtm(:) = Mtr1_mtm(:) + s1_vtw_ho_all(:,ictaper) &
        * conjg(s2_vtw_ho_all(:,ictaper))
    Mtr2_mtm(:) = Mtr2_mtm(:) + s2_vtw_ho_all(:,ictaper) &
        * conjg(s1_vtw_ho_all(:,ictaper))

    Mtr3_mtm(:) = Mtr3_mtm(:) + s1_dtw_ho_all(:,ictaper) &
        * conjg(s2_dtw_ho_all(:,ictaper))
    Mtr4_mtm(:) = Mtr4_mtm(:) + s2_dtw_ho_all(:,ictaper) &
        * conjg(s1_dtw_ho_all(:,ictaper))
    Mtr5_mtm(:) = Mtr5_mtm(:) + s2_dtw_ho_all(:,ictaper) &
        * conjg(s2_dtw_ho_all(:,ictaper))
    enddo ! ictaper

    ! compute p_j, q_j, P_j, Q_j and adjoint source fp, fq
    fp1(1:npts) = 0.
    fq1(1:npts) = 0.
    fp2(1:npts) = 0.
    fq2(1:npts) = 0.

    do ictaper = 1,ntaper

    ! compute p_j(w) and q_j(w)
    p1_adj(:) = cmplx(0.,0.)
    pw1_adj(:) = cmplx(0.,0.)
    p2_adj(:) = cmplx(0.,0.)
    pw2_adj(:) = cmplx(0.,0.)
    q1_adj(:) = cmplx(0.,0.)
    qw1_adj(:) = cmplx(0.,0.)
    q2_adj(:) = cmplx(0.,0.)
    qw2_adj(:) = cmplx(0.,0.)

    fp1_adj(:) = cmplx(0.,0.)
    fp2_adj(:) = cmplx(0.,0.)
    fq1_adj(:) = cmplx(0.,0.)
    fq2_adj(:) = cmplx(0.,0.)

    do i = i_fstart,i_fend
    !p1_adj(i)  =  -0.5* conjg(s2_vtw_ho_all(i,ictaper)) / Mtr1_mtm(i)
    pw1_adj(i) = - 0.5 * s2_vtw_ho_all(i,ictaper) / Mtr2_mtm(i)

    !p2_adj(i)  =  0.5 * conjg(s1_vtw_ho_all(i,ictaper)) / Mtr2_mtm(i)
    pw2_adj(i) =  0.5 * s1_vtw_ho_all(i,ictaper) / Mtr1_mtm(i)

    q1_adj(i)  =  0.5 * conjg(s2_dtw_ho_all(i,ictaper)) / Mtr3_mtm(i)
    !   qw1_adj(i) =  0.5 * s2_dtw_ho_all(i,ictaper) / Mtr4_mtm(i)
    !
    q2_adj(i)  =  0.5 * conjg(s1_dtw_ho_all(i,ictaper)) / Mtr4_mtm(i) &
        - conjg(s2_dtw_ho_all(i,ictaper)) / Mtr5_mtm(i)
    !   qw2_adj(i) =  0.5 * s1_dtw_ho_all(i,ictaper) / Mtr3_mtm(i) &
    !                 - s2_dtw_ho_all(i,ictaper) / Mtr5_mtm(i)
    enddo

    ! compute P_j(w) and Q_j(w)
    ! NOTE: the MT measurement is incorporated here
    !             also note that wp_taper and wq_taper can contain
    !             uncertainty estimations
    ! adds misfit measurement dtau, dlnA
    fp1_adj(i_fstart: i_fend) &
        ! = (conjg(p1_adj(i_fstart: i_fend)) +
    =  2.0*(pw1_adj(i_fstart:i_fend)) &
        * dcmplx(ddtau_w(i_fstart: i_fend),0.) &
        * dcmplx(wp_taper(i_fstart: i_fend),0.)
    fp2_adj(i_fstart: i_fend) &
        ! = (conjg(p2_adj(i_fstart: i_fend)) +
    = 2.0*(pw2_adj(i_fstart: i_fend)) &
        * dcmplx(ddtau_w(i_fstart: i_fend),0.) &
        * dcmplx(wp_taper(i_fstart: i_fend),0.)
    fq1_adj(i_fstart: i_fend) &
        = (q1_adj(i_fstart: i_fend) + conjg(q1_adj(i_fstart: i_fend)))&
        * dcmplx(ddlnA_w(i_fstart: i_fend),0.) &
        * dcmplx(wq_taper(i_fstart: i_fend),0.)
    fq2_adj(i_fstart: i_fend) &
        = (q2_adj(i_fstart: i_fend) + conjg(q2_adj(i_fstart: i_fend)))&
        * dcmplx(ddlnA_w(i_fstart: i_fend),0.) &
        * dcmplx(wq_taper(i_fstart: i_fend),0.)

    ! IFFT into the time domain
    call fftinv(LNPT,fp1_adj,dble(REVERSE_FFT),dble(deltat),fp1_adj_t)
    call fftinv(LNPT,fp2_adj,dble(REVERSE_FFT),dble(deltat),fp2_adj_t)
    call fftinv(LNPT,fq1_adj,dble(REVERSE_FFT),dble(deltat),fq1_adj_t)
    call fftinv(LNPT,fq2_adj,dble(REVERSE_FFT),dble(deltat),fq2_adj_t)

    ! create adjoint source
    ! applies taper to time signal
    fp1(1:nlen) = fp1(1:nlen) + real(tas(1:nlen,ictaper) * &
        fp1_adj_t(1:nlen),CUSTOM_REAL)
    fp2(1:nlen) = fp2(1:nlen) + real(tas(1:nlen,ictaper) * &
        fp2_adj_t(1:nlen),CUSTOM_REAL)
    fq1(1:nlen) = fq1(1:nlen) + real(tas(1:nlen,ictaper) * &
        fq1_adj_t(1:nlen),CUSTOM_REAL)
    fq2(1:nlen) = fq2(1:nlen) + real(tas(1:nlen,ictaper) * &
        fq2_adj_t(1:nlen),CUSTOM_REAL)
    enddo  ! ictaper
    deallocate(tas)
    deallocate(s1_dtw_ho_all)
    deallocate(s1_vtw_ho_all)
    deallocate(s2_dtw_ho_all)
    deallocate(s2_vtw_ho_all)

end subroutine mtm_DD_adj

!==============================================================================
!        subroutines used in mtm_measure() and mtm_adj()
!==============================================================================

!-----------------------------------------------------------------------
subroutine CC_similarity(d1,d2,npts,&
        i_tstart1,i_tend1,i_tstart2,i_tend2,&
        window_type,&
        cc_max)
    !! measure the similarity of two waveforms based on cross-correlatiobs

    use constants
    implicit none

    ! inputs & outputs
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: d1,d2
    integer, intent(in) :: i_tstart1,i_tend1,i_tstart2,i_tend2
    integer, intent(in) :: npts,window_type
    real(kind=CUSTOM_REAL), intent(out) :: cc_max

    ! window
    integer :: nlen1,nlen2,nlen
    real(kind=CUSTOM_REAL), dimension(npts) :: d1_tw,d2_tw
    ! cc
    integer :: ishift
    real(kind=CUSTOM_REAL) :: dlnA

    !! initialization
    cc_max=0.0

    !! window
    call cc_window(d1,npts,window_type,i_tstart1,i_tend1,0,0.d0,nlen1,d1_tw)
    call cc_window(d2,npts,window_type,i_tstart2,i_tend2,0,0.d0,nlen2,d2_tw)
    if(nlen1<1 .or. nlen1>npts) print*,'check nlen1 ',nlen1
    if(nlen2<1 .or. nlen2>npts) print*,'check nlen2 ',nlen2
    nlen = max(nlen1,nlen2)
    !! cc
    call xcorr_calc(d1_tw,d2_tw,npts,1,nlen,ishift,dlnA,cc_max) !T(d1-d2)

end subroutine cc_similarity

! ---------------------------------------------------------------------------
subroutine frequency_limit (dat,nlen,deltat,i_fstart,i_fend)
    use constants
    implicit none

    integer, intent(in) :: nlen
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: dat
    real(kind=CUSTOM_REAL), intent(in) ::  deltat
    integer, intent(out) :: i_fstart,i_fend

    integer :: i,fnum
    real(kind=CUSTOM_REAL) :: df,ampmax, wtr_use
    integer :: i_ampmax, i_fstart_stop, i_fend_stop
    complex (SIZE_DOUBLE), dimension(NPT) :: dat_dtwo

    !! find reliable frequency limits
    ! fft of untapered inputs
    df = 1./(NPT*deltat)
    dat_dtwo = cmplx(0.,0.)
    dat_dtwo(1:nlen) = cmplx(real(dat(1:nlen),SIZE_DOUBLE),0.0_SIZE_DOUBLE)
    call fft(LNPT,dat_dtwo,dble(FORWARD_FFT),dble(deltat))

    fnum = NPT/2+1
    ! water level based untapered dat1
    ampmax = 0.
    i_ampmax = 1
    do i = 1, fnum   ! loop over frequencies
    if( abs(dat_dtwo(i)) > ampmax) then
        ampmax =  abs(dat_dtwo(i))
        i_ampmax = i
    endif
    enddo
    wtr_use = ampmax * WTR

    ! i_fend
    i_fend = fnum
    i_fend_stop = 0
    do i = 1,fnum
    if( abs(dat_dtwo(i)) <= abs(wtr_use) .and. i_fend_stop==0 .and. i > i_ampmax ) then
        i_fend_stop = 1
        i_fend = i
    endif
    if( abs(dat_dtwo(i)) >= 10.*abs(wtr_use) .and. i_fend_stop==1 .and. i > i_ampmax) then
        i_fend_stop = 0
        i_fend = i
    endif
    enddo

    ! i_fstart
    i_fstart = 1
    i_fstart_stop = 0
    do i = fnum,1,-1
    if( abs(dat_dtwo(i)) <= abs(wtr_use) .and. i_fstart_stop==0 .and. i < i_ampmax ) then
        i_fstart_stop = 1
        i_fstart = i
    endif
    if( abs(dat_dtwo(i)) >= 10.*abs(wtr_use) .and. i_fstart_stop==1 .and. i < i_ampmax) then
        i_fstart_stop = 0
        i_fstart = i
    endif
    enddo

    if( DISPLAY_DETAILS) then
        open(1,file=trim(output_dir)//'/spectrum',status='unknown')
        do  i =  1,2*i_fend
        write(1,'(f15.5,2e15.5)') (i-1)*df,abs(dat_dtwo(i)),wtr_use
        enddo
        close(1)
    endif

end subroutine frequency_limit
! ---------------------------------------------------------------------------
subroutine cc_error(data_dtw,syn_dtw,npts,deltat,nlen,ishift,dlnA,sigma_dt,sigma_dlnA)
    ! CHT: Estimate the uncertainty in the CC measurement
    !      based on the integrated waveform difference between the data
    !      and the reconstructed synthetics.
    ! NOTE: We implement the exact equations that are in the Latex notes.

    use constants
    implicit none
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: data_dtw, syn_dtw
    integer, intent(in) :: npts,nlen, ishift
    real(kind=CUSTOM_REAL), intent(in) :: deltat, dlnA
    real(kind=CUSTOM_REAL), intent(inout) :: sigma_dt, sigma_dlnA

    real(kind=CUSTOM_REAL), dimension(npts) :: syn_dtw_cc, syn_vtw_cc
    real(kind=CUSTOM_REAL) :: sigma_dt_top, sigma_dlnA_top, sigma_dt_bot, sigma_dlnA_bot
    integer :: i,j

    ! make corrections to syn based on cc measurement
    syn_dtw_cc(:)=0.0
    do i=1,nlen
    j=i-ishift
    if(j>=1 .and. j <=nlen) then
        syn_dtw_cc(i) = syn_dtw(j) * exp(dlnA)
    endif
    enddo

    syn_vtw_cc(:)=0.0
    ! compute cc-corrected synthetic velocity
    call compute_vel(syn_dtw_cc,npts,deltat,nlen,syn_vtw_cc)

    ! estimated uncertainty in cross-correlation travltime and amplitude
    sigma_dt_top   = sum( (data_dtw(1:nlen) - syn_dtw_cc(1:nlen) )**2 )
    sigma_dt_bot   = sum( syn_vtw_cc(1:nlen)**2 )
    sigma_dlnA_top = sigma_dt_top
    sigma_dlnA_bot = sum( (syn_dtw_cc(1:nlen))**2 )/(dlnA * dlnA)
    sigma_dt       = sqrt( sigma_dt_top / sigma_dt_bot )
    sigma_dlnA     = sqrt( sigma_dlnA_top / sigma_dlnA_bot )

    ! make sure that the uncertainty estimates are not below the water level;
    ! otherwise, the adjoint sources will blow up unreasonably
    if( sigma_dt < DT_SIGMA_MIN) sigma_dt = DT_SIGMA_MIN
    if( sigma_dlnA < DLNA_SIGMA_MIN) sigma_dlnA = DLNA_SIGMA_MIN

    if(DISPLAY_DETAILS) then
        print*
        print*,'error estimation based on cc'
        print*,'sigma_dt top = ',sigma_dt_top, 'bot =', sigma_dt_bot, &
            ' ratio =', sqrt( sigma_dt_top / sigma_dt_bot )
        print*,'sigma_dlnA top = ',sigma_dlnA_top, 'bot =', sigma_dlnA_bot,&
            ' ratio =', sqrt( sigma_dlnA_top / sigma_dlnA_bot )
        print *, 'estimate sigma_dt   : ', sigma_dt
        print *, 'estimate sigma_dlnA : ', sigma_dlnA
        open(1,file=trim(output_dir)//'/cc_error.dat',status='unknown')
        do i = 1,nlen
        write(1,'(I5,4e15.5)') i, data_dtw(i), syn_dtw(i), syn_dtw_cc(i), syn_vtw_cc(i)
        enddo
        close(1)
    endif

end subroutine cc_error
! ---------------------------------------------------------------------------
subroutine write_phase(phase_func, wvec, i_right, tshift, phi_wt, dtau_wt)
    use constants
    implicit none

    ! The transfer function maps the data2 to the CC-corrected data1;

    complex (CUSTOM_REAL), dimension(*), intent(in) :: phase_func
    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: wvec
    real(kind=CUSTOM_REAL), intent(in) ::  tshift
    integer, intent(in) ::  i_right
    real(kind=CUSTOM_REAL), dimension(*), intent(out) :: phi_wt, dtau_wt

    integer :: i, j
    real(kind=CUSTOM_REAL) :: smth, smth1, smth2 ! f0

    phi_wt(1:NPT) = 0.

    ! loop to calculate phase and amplitude
    do i = 1, i_right
    phi_wt(i) = atan2( aimag(phase_func(i)) , real(phase_func(i)) )
    enddo

    ! NOTE: the CC measurements dT (tshift) and dlnA are BOTH included
    dtau_wt(1) = tshift
    do i = 1, i_right

    if (i > 1 .and. i < i_right) then
        ! check the smoothness (2nd-order derivative) by 2*pi changes
        smth  =  phi_wt(i+1) + phi_wt(i-1) - 2.0 * phi_wt(i)
        smth1 = (phi_wt(i+1) + TWOPI) + phi_wt(i-1) - 2.0 * phi_wt(i)
        smth2 = (phi_wt(i+1) - TWOPI) + phi_wt(i-1) - 2.0 * phi_wt(i)
        if(abs(smth1).lt.abs(smth).and.abs(smth1).lt.abs(smth2).and. abs(phi_wt(i) - phi_wt(i+1)) > PHASE_STEP)then
            if (DISPLAY_DETAILS) print *, 'phase correction : 2 pi', i, phi_wt(i) - phi_wt(i+1)
            do j = i+1, i_right
            phi_wt(j) = phi_wt(j) + TWOPI
            enddo
        endif
        if(abs(smth2).lt.abs(smth).and.abs(smth2).lt.abs(smth1).and. abs(phi_wt(i) - phi_wt(i+1)) > PHASE_STEP)then
            if (DISPLAY_DETAILS) print *, 'phase correction : - 2 pi', i, phi_wt(i) - phi_wt(i+1)
            do j = i+1, i_right
            phi_wt(j) = phi_wt(j) - TWOPI
            enddo
        endif
    endif

    ! add the CC measurements to the transfer function
    if (i > 1) dtau_wt(i) = (0.5/wvec(i)) * phi_wt(i) + tshift

    enddo

end subroutine write_phase
! ---------------------------------------------------------------------------
subroutine write_trans(trans, wvec,i_left, i_right, tshift, dlnA,phi_wt, abs_wt, dtau_wt, dlnA_wt)
    ! The transfer function maps the synthetics to the CC-deconstructed data;
    ! the CC measurements then need to be applied to match the original data.

    use constants
    implicit none
    complex (CUSTOM_REAL), dimension(*),intent(in) :: trans
    real(kind=CUSTOM_REAL), dimension(*),intent(in) :: wvec
    real(kind=CUSTOM_REAL), intent(in) :: tshift, dlnA
    integer, intent(in) :: i_left, i_right
    real(kind=CUSTOM_REAL), dimension(*), intent(out) :: phi_wt, abs_wt,dtau_wt,dlnA_wt

    integer :: i

    ! initialization
    phi_wt(1:NPT) = 0.d0
    abs_wt(1:NPT) = 0.d0
    dtau_wt(1:NPT) = 0.d0
    dlnA_wt(1:NPT) = 0.d0

    ! loop to calculate phase and amplitude
    do i = i_left, i_right
    phi_wt(i) = atan2( aimag(trans(i)) , real(trans(i)) )
    abs_wt(i) = abs(trans(i))
    enddo

    ! check the smoothenss of phi
    call smooth_phase(phi_wt,i_left,i_right)

    ! NOTE: the CC measurements dT (tshift) and dlnA are BOTH included
    ! print*,' add the CC measurements to the transfer function'
    if (i_left > 1) dtau_wt(i_left:i_right) = (-1./wvec(i_left:i_right)) * phi_wt(i_left:i_right) + tshift
    if (i_left == 1) dtau_wt(i_left+1:i_right) = (-1./wvec(i_left+1:i_right)) * phi_wt(i_left+1:i_right) + tshift
    dtau_wt(1) = tshift
    dlnA_wt(i_left:i_right) = log(abs_wt(i_left:i_right)) + dlnA

end subroutine write_trans
! --------------------------------------------------------------------
subroutine smooth_phase(phi_wt,i_left,i_right)
    ! check the smoothness of phi
    use constants
    implicit none

    real(kind=CUSTOM_REAL), dimension(*),intent(out) :: phi_wt
    integer, intent(in) :: i_left, i_right

    integer :: i, j
    real(kind=CUSTOM_REAL) :: smth, smth1, smth2 ! f0

    do i = i_left, i_right
    if (i > i_left .and. i < i_right) then
        !  print*,' check the smoothness (2nd-order derivative) by 2*pi changes'
        smth  =  phi_wt(i+1) + phi_wt(i-1) - 2.0 * phi_wt(i)
        smth1 = (phi_wt(i+1) + TWOPI) + phi_wt(i-1) - 2.0 * phi_wt(i)
        smth2 = (phi_wt(i+1) - TWOPI) + phi_wt(i-1) - 2.0 * phi_wt(i)
        if(abs(smth1).lt.abs(smth).and.abs(smth1).lt.abs(smth2).and. abs(phi_wt(i) - phi_wt(i+1)) > PHASE_STEP)then
            if (DISPLAY_DETAILS) print *, 'phase correction : 2pi',i,phi_wt(i) - phi_wt(i+1)
            do j = i+1, i_right
            phi_wt(j) = phi_wt(j) + TWOPI
            enddo
        endif

        if(abs(smth2).lt.abs(smth).and.abs(smth2).lt.abs(smth1).and. abs(phi_wt(i) - phi_wt(i+1)) > PHASE_STEP)then
            if (DISPLAY_DETAILS) print *, 'phase correction : - 2 pi ', i, phi_wt(i) - phi_wt(i+1)
            do j = i+1, i_right
            phi_wt(j) = phi_wt(j) - TWOPI
            enddo
        endif

    endif
    enddo

end subroutine smooth_phase
! --------------------------------------------------------------------
subroutine unwrap_phase(phi_wt,i_left,i_right)
    ! unwrap phi
    use constants
    implicit none

    real(kind=CUSTOM_REAL), dimension(*),intent(out) :: phi_wt
    integer, intent(in) :: i_left, i_right

    integer :: i
    real(kind=CUSTOM_REAL) :: dp, dps, dp_corr

    do i = i_left, i_right
    dp_corr = 0.0
    if (i > i_left .and. i <= i_right) then
        !% Incremental phase variations
        dp = phi_wt(i)-phi_wt(i-1)
        !% Equivalent phase variations in [-pi,pi)
        dps = mod(dp+PI,TWOPI)-PI
        !% Preserve variation sign for pi vs. -pi
        if(dps==-PI .and. dp>0.0) dps=PI
        !% Incremental phase corrections
        dp_corr = dps-dp
        !% Ignore correction when incr. variation is < CUTOFF
        if(abs(dp)<CUTOFF) dp_corr=0.0
        !% Integrate corrections and add to P to produce smoothed phase values
        phi_wt(i:i_right) = phi_wt(i:i_right) + dp_corr
    endif
    enddo

end subroutine unwrap_phase
! --------------------------------------------------------------------
subroutine time_costaper(dat,syn,npts,istart,iend)
    use constants
    implicit none

    ! input parameters
    integer, intent(in) :: istart,iend,npts
    real(kind=CUSTOM_REAL), dimension(*), intent(out) :: dat,syn

    integer :: B1,B2,B3,B4,nlen_taper
    integer i, n_width_left,n_width_right,nlen
    real(kind=CUSTOM_REAL) :: fac

    ! set the number of points corresponding to the taper width
    if(istart<iend) then
        nlen_taper=floor((iend-istart)/4.0)
        B1=max(istart-nlen_taper,1)
        B4=min(iend+nlen_taper,npts)
        B2=istart
        B3=iend
    endif
    n_width_left=B2-B1
    n_width_right=B4-B3

    ! left side
    nlen=2*n_width_left
    do i = 1,n_width_left
    !fac = 1.                                         ! boxcar window
    !fac = 1 - sfac1*((i-1) - real(nlen)/2.)**2       ! welch window
    fac = 1. - cos(PI*(i-1)/(nlen-1))**ipwr_t
    dat(B1+i-1)=dat(B1-1+i)*fac
    syn(B1+i-1)=syn(B1-1+i)*fac
    enddo

    ! right side
    nlen=2*n_width_right
    do i = 1,n_width_right
    !fac = 1.                                         ! boxcar window
    !fac = 1 - sfac1*((i-1) - real(nlen)/2.)**2       ! welch window
    fac = 1. - cos(PI*(i-1)/(nlen-1))**ipwr_t
    dat(B4-i+1)=dat(B4-i+1)*fac
    syn(B4-i+1)=syn(B4-i+1)*fac
    enddo

    ! beyond the window is zero
    dat(1:B1)=0.0;
    dat(B4:npts)=0.0
    syn(1:B1)=0.0;
    syn(B4:npts)=0.0

end subroutine time_costaper
! ---------------------------------------------------------------------------
subroutine cc_window(dat,npts,window_type,istart,iend,ishift,dlnA,nlen,dat_win)
    ! delay by ishift and scale by exp(dlnA)
    use constants
    implicit none

    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: dat
    real(kind=CUSTOM_REAL), intent(in) :: dlnA
    integer, intent(in) :: npts, istart,iend, ishift,window_type
    integer, intent(out) :: nlen
    real(kind=CUSTOM_REAL), dimension(*), intent(out) :: dat_win

    integer ::  i, j
    real(kind=CUSTOM_REAL) :: sfac1,  fac

    !print*
    !print*,'window with corrections: isfhit = ', ishift, ' dlnA = ',dlnA,&
    !'window type -- ',window_type

    nlen = iend - istart+1

    ! some constants
    sfac1 = (2./real(nlen))**2   ! for Welch window

    ! initialization
    dat_win(1:npts) = 0.d0

    do i = 1, nlen
    if(window_type ==2) then
        fac = 1 - sfac1*((i-1) - real(nlen)/2.)**2
    elseif(window_type ==3) then
        fac = 1. - cos(PI*(i-1)/(nlen-1))**ipwr_t
    elseif(window_type ==4) then
        fac = 0.5 - 0.5*cos(TWOPI*(i-1)/(nlen-1))
    else
        fac = 1. ! boxcar window
    endif

    ! index
    j = i + (istart-1) - ishift

    if(j>=1 .and. j <=npts) then
        dat_win(i) = dat(j) * exp(dlnA) * fac ! shift and scale
    endif

    enddo

end subroutine cc_window
! ---------------------------------------------------------------------------
subroutine cc_window_inverse(dat_win,npts,window_type,istart,iend,ishift,dlnA,dat)
    use constants
    implicit none

    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: dat_win
    real(kind=CUSTOM_REAL), intent(in) :: dlnA
    integer, intent(in) :: npts, istart,iend, ishift,window_type
    real(kind=CUSTOM_REAL), dimension(*), intent(out) :: dat

    integer ::  i, j, nlen
    real :: sfac1, fac

    !  print*
    !  print*,'window with corrections: isfhit = ', ishift, ' dlnA = ',dlnA

    nlen = iend - istart

    ! some constants
    sfac1 = (2./real(nlen))**2   ! for Welch window

    ! initialization
    dat(1:npts) = 0.d0

    do i = 1, nlen
    if(window_type ==2) then
        fac = 1 - sfac1*((i-1) - real(nlen)/2.)**2
    else if(window_type ==3) then
        fac = 1. - cos(PI*(i-1)/(nlen-1))**ipwr_t
    else if(window_type ==4) then
        fac = 0.5 - 0.5*cos(TWOPI*(i-1)/(nlen-1))
    else
        fac = 1. ! boxcar window
    endif

    j = i + (istart-1) - ishift

    ! window dat2 using exact window info
    if(j>=1 .and. j <=npts) then
        dat(j) = dat_win(i) * exp(-dlnA) * fac
    endif
    enddo
end subroutine cc_window_inverse
! ---------------------------------------------------------------------------
subroutine compute_cc(dat1, dat2, nlen, dt, ishift, tshift, dlnA, cc_max)
    use constants
    implicit none

    ! time shift MEASUREMENT between data1 and data2
    ! CHT: modified the subroutine to resemble the one used in FLEXWIN
    ! tshift = T(dat1) - T(dat2)
    ! dlnA = 0.5 * log(A(dat1)**2 / A(dat2)**2)

    real(kind=CUSTOM_REAL), dimension(*), intent(in) :: dat1, dat2
    integer, intent(in) :: nlen
    real(kind=CUSTOM_REAL), intent(in) :: dt
    real(kind=CUSTOM_REAL), intent(out) :: tshift, dlnA, cc_max
    integer, intent(out) :: ishift

    real(kind=CUSTOM_REAL) :: cc, norm_s, norm ! cr_shift
    integer i1, i2, i, j, i_left, i_right, id_left, id_right
    ! initialise shift and cross correlation to zero
    ishift = 0
    cc_max = 0.0

    ! index of window limits
    i1 = 1
    i2 = nlen

    ! length of window (number of points, including ends)
    !nlen = i2 - i1 + 1

    ! power of synthetic signal in window
    norm_s = sqrt(sum(dat2(i1:i2)*dat2(i1:i2)))

    ! left and right limits of index (time) shift search
    ! NOTE: This looks OUTSIDE the time window of interest to compute TSHIFT and
    ! CC.
    !       How far to look outside, in theory, should be another parameter.
    !       However, it does not matter as much if the data and synthetics are
    !          zeroed outside the windows.
    i_left = - nlen
    i_right = nlen

    ! i is the index to shift to be applied to DATA (data)
    do i = i_left, i_right

    ! normalization factor varies as you take different windows of data
    id_left = max(1,i1+i)      ! left index for data window
    id_right = min(nlen,i2+i)  ! right index for data window
    norm = norm_s * sqrt(sum(dat1(id_left:id_right)*(dat1(id_left:id_right))))

    ! cc as a function of i
    cc = 0.
    do j = i1, i2   ! loop over full window length
    if((j+i).ge.1 .and. (j+i).le.nlen) cc = cc + dat2(j)*dat1(j+i)  ! d ishifted by i
    enddo
    cc = cc/norm

    if (cc > cc_max) then
        ! CHT: do not allow time shifts larger than the specified input range
        ! This is an important criterion, since it may pick TSHIFT_MIN or
        ! TSHIFT_MAX
        ! if cc_max within the interval occurs on the boundary.
        !         if( (i*dt >= TSHIFT_MIN).and.(i*dt <= TSHIFT_MAX) ) then
        cc_max = cc
        ishift = i
        !         endif
    endif

    enddo
    tshift = ishift*dt
    ! The previously used expression for dlnA of Dahlen and Baig (2002),
    ! is a first-order perturbation of ln(A1/A2) = (A1-A2)/A2 .
    ! The new expression is better suited to getting Gaussian-distributed
    ! values between -1 and 1, with dlnA = 0 indicating perfect fit, as before.
    dlnA = 0.5 * log( sum(dat1(i1:i2) * dat1(i1:i2)) / sum(dat2(i1:i2) * dat2(i1:i2)) )

end subroutine compute_cc
!-----------------------------------------------------------------------------
subroutine compute_vel(syn,npts,deltat,nlen,syn_vel)
    use constants
    implicit none
    real(kind=CUSTOM_REAL), dimension(*),intent(in) :: syn
    real(kind=CUSTOM_REAL), intent(in) :: deltat
    integer, intent(in) :: npts, nlen
    real(kind=CUSTOM_REAL), dimension(*),intent(out) :: syn_vel

    ! index
    integer :: itime

    ! initialization
    syn_vel(1:npts)=0.0

    do itime=2,nlen-1
    syn_vel(itime)=(syn(itime+1)-syn(itime-1))/(2*deltat)
    enddo
    ! boundaries
    syn_vel(1)=(syn(2)-syn(1))/deltat
    syn_vel(nlen)=(syn(nlen)-syn(nlen-1))/deltat

end subroutine compute_vel
!-----------------------------------------------------------------------------


module derivsmod
!$$$   module documentation block
!                .      .    .                                       .
! module:    derivsmod
!
! abstract:  This module defines and holds bundles that contain
!            guess derivative fields.
!
! program history log:
!   2013-10-19 Todling - Initial code.
!   2014-06-18 Carley - add lgues and dlcbasdlog
!   2015-07-10 Pondeca - add cldchgues and dcldchdlog
!   2016-05-10 Thomas - remove references to cwgues0
!   2019-05-08 mtong - replace set_ with init_anadv
!   2019-05-08 eliu - recover logic (drv_set) to indicate the derivative
!                     vars are allocated and defined
!
! public subroutines:
!  drv_initialized         - initialize name of fields to calc derivs for
!  create_ges_derivatives  - initialize bundles holding derivatives
!  destroy_ges_derivatives - finalize bundles holding derivatives

! public variables:
!  gsi_xderivative_bundle  - bundle holding longitudinal derivatives
!  gsi_yderivative_bundle  - bundle holding latitudinal derivatives
!  dvars2d, dvars3d        - names of 2d/3d derivatives
!  dsrcs2d, dsrcs3d        - names of where original fields reside
!  drv_initialized         - flag indicating initialization status
!  drv_set                 - flag indicating the variables are allocated and defined
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP; SGI Origin 2000; Compaq/HP
!
!$$$ end documentation block

use kinds, only: i_kind, r_kind
use mpimod, only: mype
use gridmod, only: lat2,lon2,nsig
use constants, only: zero,max_varname_length
use state_vectors, only: svars2d,svars3d
use GSI_BundleMod, only : GSI_BundleCreate
use GSI_BundleMod, only : GSI_Bundle
use GSI_BundleMod, only : GSI_BundleGetPointer
use GSI_BundleMod, only : GSI_BundleDestroy
use GSI_BundleMod, only : GSI_BundleUnset

use GSI_BundleMod, only : GSI_Grid
use GSI_BundleMod, only : GSI_GridCreate

use GSI_MetGuess_Mod, only: gsi_metguess_bundle
use GSI_ChemGuess_Mod, only: gsi_chemguess_bundle

use mpeu_util, only: getindex
implicit none
save
private

public :: drv_initialized
public :: drv_set
public :: create_ges_derivatives
public :: destroy_ges_derivatives

public :: gsi_xderivative_bundle
public :: gsi_yderivative_bundle
public :: dvars2d, dvars3d
public :: dsrcs2d, dsrcs3d
public :: cwgues,cfgues
public :: ggues,vgues,pgues,lgues,dvisdlog,dlcbasdlog
public :: w10mgues,howvgues,cldchgues,dcldchdlog
public :: qsatg,qgues,dqdt,dqdrh,dqdp
public :: init_anadv

logical :: drv_initialized = .false.
logical :: drv_set = .false.

type(gsi_bundle),pointer :: gsi_xderivative_bundle(:)
type(gsi_bundle),pointer :: gsi_yderivative_bundle(:)
character(len=max_varname_length),allocatable,dimension(:):: dvars2d, dvars3d
character(len=max_varname_length),allocatable,dimension(:):: dsrcs2d, dsrcs3d

real(r_kind),allocatable,dimension(:,:,:):: qsatg,qgues,dqdt,dqdrh,dqdp
real(r_kind),allocatable,dimension(:,:):: ggues,vgues,pgues,lgues,dvisdlog,dlcbasdlog
real(r_kind),allocatable,dimension(:,:):: w10mgues,howvgues,cldchgues,dcldchdlog
real(r_kind),target,allocatable,dimension(:,:,:):: cwgues,cfgues

! below this point: declare vars not to be made public

character(len=*),parameter:: myname='derivsmod'
integer(i_kind),allocatable,dimension(:):: levels
contains

subroutine init_anadv
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:	 define derivatives
!   prgmmr:	 todling
!      org:	 NASA/GSFC, Global Modeling and Assimilation Office, 610.1
!     date:	 2013-10-19
!
! abstract: - set which derivatives to calculate from guess fields
!
! program history log:
!   2013-09-27  todling  - initial code
!   2014-02-03  todling  - negative levels mean rank-3 array
!   2019-05-08  mtong    - replace set_ with init_anadv
!   2019-05-08  eliu     - recover logic (drv_set) to indicate the derivative
!                          vars are allocated and defined
!
!   input argument list: see Fortran 90 style document below
!
!   output argument list: see Fortran 90 style document below
!
! attributes:
!   language: Fortran 90 and/or above
!   machine:
!
!$$$  end subprogram documentation block
use mpeu_util, only: gettablesize
use mpeu_util, only: gettable
use mpeu_util, only: getindex
implicit none

character(len=*),parameter:: rcname='anavinfo'

character(len=*),parameter::myname_=myname//'*set_'
character(len=*),parameter:: tbname='state_derivatives::'
integer(i_kind) luin,ii,nrows,ntot,ipnt,istatus
integer(i_kind) i2d,i3d,n2d,n3d,irank
integer(i_kind),allocatable,dimension(:)::nlevs
character(len=256),allocatable,dimension(:):: utable
character(len=max_varname_length),allocatable,dimension(:):: vars
character(len=max_varname_length),allocatable,dimension(:):: sources
logical matched

if(drv_set) return

open(newunit=luin,file=trim(rcname),form='formatted')

! Scan file for desired table first
! and get size of table
call gettablesize(tbname,luin,ntot,nrows)
if(nrows==0) then
   if(luin/=5) close(luin)
   return
endif

! Get contents of table
allocate(utable(nrows))
call gettable(tbname,luin,ntot,nrows,utable)

! release file unit
if(luin/=5) close(luin)

! allocate space for entries from table
allocate(vars(nrows),nlevs(nrows),sources(nrows))

! Retrieve each token of interest from table and define
! variables participating in state vector
n2d=0; n3d=0
do ii=1,nrows
   read(utable(ii),*) vars(ii),&  ! variable name
                      nlevs(ii),& ! number of levels
                      sources(ii) ! source
   if (nlevs(ii)==1) then
      n2d=n2d+1
   else
      n3d=n3d+1
   endif
enddo

deallocate(utable)

allocate(dvars2d(n2d),dvars3d(n3d),&
         dsrcs2d(n2d),dsrcs3d(n3d),levels(n3d))

! loop over variables and identify them by comparison
i2d=0; i3d=0
do ii=1,nrows
   matched=.false.
   if(trim(sources(ii))=='met_guess') then
      if(associated(gsi_metguess_bundle)) then
         call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(vars(ii)),ipnt,istatus,irank=irank);
         if (ipnt>0) then
            if(irank==2) then
               i2d=i2d+1
               dvars2d(i2d)=trim(vars(ii))
               dsrcs2d(i2d)=trim(sources(ii))
               matched=.true.
            endif
            if(irank==3) then
               i3d=i3d+1
               dvars3d(i3d)=trim(vars(ii))
               dsrcs3d(i3d)=trim(sources(ii))
               levels(i3d) =abs(nlevs(ii))
               matched=.true.
            endif
         endif
      endif
   endif
   if(trim(sources(ii))=='chem_guess') then
      if(associated(gsi_chemguess_bundle)) then
         call gsi_bundlegetpointer(gsi_chemguess_bundle(1),trim(vars(ii)),ipnt,istatus,irank=irank);
         if (ipnt>0) then
            if(irank==2) then
               i2d=i2d+1
               dvars2d(i2d)=trim(vars(ii))
               dsrcs2d(i2d)=trim(sources(ii))
               matched=.true.
            endif
            if(irank==3) then
               i3d=i3d+1
               dvars3d(i3d)=trim(vars(ii))
               dsrcs3d(i3d)=trim(sources(ii))
               levels(i3d) =abs(nlevs(ii))
               matched=.true.
            endif
         endif
      endif
   endif
   ! now care for variables not in guess (usually, derived tendencies)
   if (.not.matched) then
      if(nlevs(ii)==1) then
         i2d=i2d+1
         dvars2d(i2d)=trim(vars(ii))
         dsrcs2d(i2d)='derived'
      else
         i3d=i3d+1
         dvars3d(i3d)=trim(vars(ii))
         dsrcs3d(i3d)='derived'
         levels(i3d) =abs(nlevs(ii))
      endif
   endif
enddo

if (mype == 0) then
    write(6,*) myname_,':  DERIVATIVE VARIABLES: '
    write(6,*) myname_,':  2D-DERV STATE VARIABLES: '
    do ii=1,n2d
       write(6,*) trim(dvars2d(ii))
    enddo
    write(6,*) myname_,':  3D-DERV STATE VARIABLES:'
    do ii=1,n3d
       write(6,*) trim(dvars3d(ii))
    enddo
end if

deallocate(vars,nlevs,sources)
drv_set=.true.

 end subroutine init_anadv

 subroutine create_ges_derivatives(switch_on_derivatives,nfldsig)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:	 create derivatives
!   prgmmr:	 todling
!      org:	 NASA/GSFC, Global Modeling and Assimilation Office, 610.1
!     date:	 2013-10-19
!
! abstract: - allocate space for bundle-kept derivatives
!
! program history log:
!   2013-09-27  todling  - initial code
!
!   input argument list: see Fortran 90 style document below
!
!   output argument list: see Fortran 90 style document below
!
! attributes:
!   language: Fortran 90 and/or above
!   machine:
!
!$$$  end subprogram documentation block
 implicit none

 logical, intent(in) :: switch_on_derivatives
 integer(i_kind),intent(in) :: nfldsig

 integer nt,ierror
 character(len=32) bname
 type(gsi_grid) :: grid

! Extra mambo-jambo
  call create_auxiliar_

  if (.not.switch_on_derivatives) return
  if (drv_initialized) return

! create derivative grid
  call GSI_GridCreate(grid,lat2,lon2,nsig)

! allocate structures
  allocate(gsi_xderivative_bundle(nfldsig))
  allocate(gsi_yderivative_bundle(nfldsig))

! Note:
!   Original code needed ps derivatives in all time slots
!   Present code creates all derivatives in all time slots
  do nt=1,nfldsig

!    create logitudinal derivative bundle
     write(bname,'(a,i3.3)') 'Lon Derivative Vector-',nt
     call GSI_BundleCreate(gsi_xderivative_bundle(nt),grid,bname,ierror, &
                           names2d=dvars2d,names3d=dvars3d,levels=levels,bundle_kind=r_kind)

!    create latidutinal derivative bundle
     write(bname,'(a,i3.3)') 'Lat Derivative Vector-',nt
     call GSI_BundleCreate(gsi_yderivative_bundle(nt),grid,bname,ierror, &
                           names2d=dvars2d,names3d=dvars3d,levels=levels,bundle_kind=r_kind)

  enddo

  drv_initialized = .true.

  if(mype==0) write(6,*) 'create_ges_derivatives: successfully complete'
  end subroutine create_ges_derivatives

  subroutine destroy_ges_derivatives
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:	 destroy derivatives
!   prgmmr:	 todling
!      org:	 NASA/GSFC, Global Modeling and Assimilation Office, 610.1
!     date:	 2013-10-19
!
! abstract: - deallocates bundles keeping derivatives
!
! program history log:
!   2013-09-27  todling  - initial code
!
!   input argument list: see Fortran 90 style document below
!
!   output argument list: see Fortran 90 style document below
!
! attributes:
!   language: Fortran 90 and/or above
!   machine:
!
!$$$  end subprogram documentation block
  use mpimod, only: mpi_comm_world
  implicit none
  integer(i_kind) nt,ierror

  if(.not.drv_initialized) return

! destroy mambo-jambo
  call destroy_auxiliar_

! destroy each instance of derivatives
  do nt=1,size(gsi_yderivative_bundle)

!    create logitudinal derivative bundle
     call GSI_BundleDestroy(gsi_yderivative_bundle(nt),ierror)
     if(ierror/=0) then
        if(mype==0) write(6,*)'warning: y-derivative not properly destroyed'
     endif

!    create latidutinal derivative bundle
     call GSI_BundleDestroy(gsi_xderivative_bundle(nt),ierror)
     if(ierror/=0) then
        if(mype==0) write(6,*)'warning: x-derivative not properly destroyed'
     endif

  enddo

! deallocate structures
  deallocate(gsi_xderivative_bundle)
  deallocate(gsi_yderivative_bundle)

! destroy derivative grid
! call GSI_GridDestroy(grid,lat2,lon2,nsig)

  deallocate(dvars2d,dvars3d,&
             dsrcs2d,dsrcs3d,levels)

  if(mype==0) write(6,*) 'destroy_ges_derivatives: successfully complete'
  end subroutine destroy_ges_derivatives

  subroutine create_auxiliar_
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    create_auxiliar_
!   prgmmr: treadon          org: np23               date:  2003-11-24
!
! abstract: allocate memory for cost function variables
!
! program history log:
!   2003-11-24  treadon
!   2004-05-18  kleist, documentation
!   2004-07-28  treadon - simplify subroutine argument list
!   2005-03-28  wu - replace mlath with mlat, modify dim of varq
!   2005-06-15  treadon - remove "use guess_grids"
!   2008-05-12  safford - rm unused uses
!   2011-02-16  zhu     - add ggues,vgues,pgues
!   2011-07-15  zhu     - add cwgues
!   2013-10-25  todling - revisit variable initialization
!   2013-11-12  lueken - revisit logic around cwgues
!   2014-02-03  todling - CV length and B-dims here (no longer in observer)
!   2014-03-19  pondeca - add w10mgues
!   2014-05-07  pondeca - add howvgues
!   2014-06-18  carley - add lgues and dlcbasdlog
!   2015-07-10  pondeca- add cldchgues and dcldchdlog
!
!   input argument list:
!    mlat
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    use constants, only: zero
    use gridmod, only: lat2,lon2,nsig
    implicit none

    integer(i_kind) i,j,k

    if (getindex(svars3d,'q')>0) then
       allocate(qsatg(lat2,lon2,nsig),&
            dqdt(lat2,lon2,nsig),dqdrh(lat2,lon2,nsig),&
            dqdp(lat2,lon2,nsig),&
            qgues(lat2,lon2,nsig))

       do k=1,nsig
          do j=1,lon2
             do i=1,lat2
                qsatg(i,j,k)=zero
                dqdt(i,j,k)=zero
                dqdrh(i,j,k)=zero
                dqdp(i,j,k)=zero
                qgues(i,j,k)=zero
             end do
          end do
       end do
    endif

    allocate(cwgues(lat2,lon2,nsig))
    do k=1,nsig
       do j=1,lon2
          do i=1,lat2
             cwgues(i,j,k)=zero
          end do
        end do
    end do

    allocate(cfgues(lat2,lon2,nsig))
    do k=1,nsig
       do j=1,lon2
          do i=1,lat2
             cfgues(i,j,k)=zero
          end do
        end do
    end do

    if (getindex(svars2d,'gust')>0) then
       allocate(ggues(lat2,lon2))
       do j=1,lon2
          do i=1,lat2
             ggues(i,j)=zero
          end do
       end do
    end if
    if (getindex(svars2d,'vis')>0) then
       allocate(vgues(lat2,lon2),dvisdlog(lat2,lon2))
       do j=1,lon2
          do i=1,lat2
             vgues(i,j)=zero
             dvisdlog(i,j)=zero
          end do
       end do
    end if
    if (getindex(svars2d,'pblh')>0) then
       allocate(pgues(lat2,lon2))
       do j=1,lon2
          do i=1,lat2
             pgues(i,j)=zero
          end do
       end do
    end if
    if (getindex(svars2d,'lcbas')>0) then
       allocate(lgues(lat2,lon2),dlcbasdlog(lat2,lon2))
       do j=1,lon2
          do i=1,lat2
             lgues(i,j)=zero
             dlcbasdlog(i,j)=zero
          end do
       end do
    end if
    if (getindex(svars2d,'wspd10m')>0) then
       allocate(w10mgues(lat2,lon2))
       do j=1,lon2
          do i=1,lat2
             w10mgues(i,j)=zero
          end do
       end do
    end if
    if (getindex(svars2d,'howv')>0) then
       allocate(howvgues(lat2,lon2))
       do j=1,lon2
          do i=1,lat2
             howvgues(i,j)=zero
          end do
       end do
    end if
    if (getindex(svars2d,'cldch')>0) then
       allocate(cldchgues(lat2,lon2),dcldchdlog(lat2,lon2))
       do j=1,lon2
          do i=1,lat2
             cldchgues(i,j)=zero
             dcldchdlog(i,j)=zero
          end do
       end do
    end if

    return
  end subroutine create_auxiliar_

  subroutine destroy_auxiliar_
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    destroy_auxliar_
!   prgmmr: treadon          org: np23               date:  2003-11-24
!
! abstract: deallocate memory from cost function variables
!
! program history log:
!   2003-11-24  treadon
!   2004-05-18  kleist, documentation
!   2011-02-16  zhu     - add ggues,vgues,pgues
!   2011-07-15  zhu     - add cwgues
!   2013-10-25  todling, revisit deallocs
!   2014-03-19  pondeca - add w10mgues
!   2014-05-07  pondeca - add howvgues
!   2014-06-18  carley - add lgues and dlcbasdlog
!   2015-07-10  pondeca- add cldchgues and dcldchdlog
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    implicit none

    if(allocated(dqdt)) deallocate(dqdt)
    if(allocated(dqdrh)) deallocate(dqdrh)
    if(allocated(dqdp)) deallocate(dqdp)
    if(allocated(qsatg)) deallocate(qsatg)
    if(allocated(qgues)) deallocate(qgues)
    if(allocated(cwgues)) deallocate(cwgues)
    if(allocated(cfgues)) deallocate(cfgues)
    if(allocated(ggues)) deallocate(ggues)
    if(allocated(vgues)) deallocate(vgues)
    if(allocated(dvisdlog)) deallocate(dvisdlog)
    if(allocated(pgues)) deallocate(pgues)
    if(allocated(lgues)) deallocate(lgues)
    if(allocated(dlcbasdlog)) deallocate(dlcbasdlog)
    if(allocated(w10mgues)) deallocate(w10mgues)
    if(allocated(howvgues)) deallocate(howvgues)
    if(allocated(cldchgues)) deallocate(cldchgues)
    if(allocated(dcldchdlog)) deallocate(dcldchdlog)

    return
  end subroutine destroy_auxiliar_
end module derivsmod


text: x <= (((a or b) and d) or a or b or c);


text:
!!
!! @HEADER
!!
!!!!**********************************************************************
!!
!!  Zoltan Toolkit for Load-balancing, Partitioning, Ordering and Coloring
!!                  Copyright 2012 Sandia Corporation
!!
!! Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
!! the U.S. Government retains certain rights in this software.
!!
!! Redistribution and use in source and binary forms, with or without
!! modification, are permitted provided that the following conditions are
!! met:
!!
!! 1. Redistributions of source code must retain the above copyright
!! notice, this list of conditions and the following disclaimer.
!!
!! 2. Redistributions in binary form must reproduce the above copyright
!! notice, this list of conditions and the following disclaimer in the
!! documentation and/or other materials provided with the distribution.
!!
!! 3. Neither the name of the Corporation nor the names of the
!! contributors may be used to endorse or promote products derived from
!! this software without specific prior written permission.
!!
!! THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
!! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
!! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
!! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
!! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
!! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
!! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
!! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
!! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
!! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!!
!! Questions? Contact Karen Devine	kddevin@sandia.gov
!!                    Erik Boman	egboman@sandia.gov
!!
!!!!**********************************************************************
!!
!! @HEADER
 !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                          //
!! File:      driver.cc                                                     //
!! Project:   Local HSFC Ordering                                           //
!! Author:    Michael Wolf                                                  //
!! Date Started:  11/02/2009                                                //
!!                                                                          //
!! Description:                                                             //
!!              File tests local HSFC ordering for simple test problem      //
!!                                                                          //
!! $Id: driver.cc 11 2009-11-10 00:15:18Z mmwolf $                          //
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program exampleRCB
   use mpi_h
   use zoltan
   use zoltanRCBex

   implicit none

   integer(Zoltan_INT) :: error
   real(Zoltan_FLOAT) :: version

!!!  numGlobObjs, numLocObjs, GIDs, xcoords, ycoords defined in zoltanRCBex module

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   interface

     subroutine readInputObjects(fname,numGlobObjs,numLocObjs,GIDs,xs,ys)
        character (len=*) :: fname
        integer :: numGlobObjs, numLocObjs
        integer, allocatable :: GIDs(:)
        real, allocatable :: xs(:), ys(:)
     end subroutine readInputObjects

   end interface
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !! Body of program
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   call MPI_Init(error)

   error = Zoltan_Initialize(version)

   call readInputObjects("mesh.txt",numGlobObjs,numLocObjs,GIDs,xcoords,ycoords)

   call partitionMeshWithRCB()

   call visualizePartition()

   deallocate(GIDs)
   deallocate(xcoords,ycoords)

   !! function in zoltanRCBex module that cleans up Zoltan data structures
   call zoltanCleanUp()

   call MPI_Finalize(error)
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

end program exampleRCB
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine readInputObjects(fname, numGlobObjs, numLocObjs,GIDs,xcoords,ycoords)
  use mpi_h
  implicit none

  character (len=*) :: fname

  integer :: numGlobObjs, numLocObjs
  integer, allocatable :: GIDs(:)
  real, allocatable :: xcoords(:), ycoords(:)

  ! Local declarations
  integer :: fnum = 2, i, currIndx
  integer :: myRank, numProcs, mpi_ierr
  integer :: tmpGID
  real :: tmpX, tmpY

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  call MPI_Comm_rank(MPI_COMM_WORLD, myrank, mpi_ierr)
  call MPI_Comm_size(MPI_COMM_WORLD, numProcs, mpi_ierr)

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! Scan data to determine # local objects
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  open(unit = fnum, file = fname)
  read(fnum,*) numGlobObjs

  numLocObjs=0

  do i = 1, numGlobObjs, 1
    read(fnum,*) tmpGID, tmpX, tmpY

    !! assumes gids start at 1, gives round robin initial distribution
    if ( MOD(tmpGID-1,numProcs) == myRank) then
      numLocObjs = numLocObjs + 1
    end if
  end do

  close(fnum)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! Allocate data for my part of mesh
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  allocate(GIDs(numLocObjs))
  allocate(xcoords(numLocObjs))
  allocate(ycoords(numLocObjs))
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! Fill data for my part of mesh
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  open(unit = fnum, file = fname)

  read(fnum,*) numGlobObjs

  currIndx = 1
  do i = 1, numGlobObjs, 1
    read(fnum,*) tmpGID, tmpX, tmpY

    !! assumes gids start at 1, gives round robin initial distribution
    if ( MOD(tmpGID-1,numProcs) == myRank) then
      GIDs(currIndx) = tmpGID
      xcoords(currIndx) = tmpX
      ycoords(currIndx) = tmpY
      currIndx = currIndx + 1
    end if
  end do

  close(fnum)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end subroutine readInputObjects
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Works for this specific 6x6 mesh
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine showSimpleMeshPartitions(myProc, numIDs, GIDs, parts)
  use mpi_h
  implicit none

  integer :: myProc, numIDs
  integer :: GIDs(*), parts(*)

  !! Local variables
  integer :: partAssign(36)
  integer :: allPartAssign(36)
  integer :: i, j, part, mpi_ierr
  integer :: partRow(6)

  data partAssign/ 36 * 0/

  do i = 1, numIDs, 1
    partAssign(GIDs(i)) = parts(i);
  end do

  call MPI_Reduce(partAssign, allPartAssign, 36, MPI_INTEGER, MPI_MAX, 0, MPI_COMM_WORLD, mpi_ierr);

  if (myProc == 0) then
    do i=30, 0, -6
      do j = 1, 6, 1
        partRow(j) = allPartAssign(i+j);
      end do

      write(*,'(I1,A,I1,A,I1,A,I1,A,I1,A,I1)') partRow(1), '-----', partRow(2), '-----', partRow(3), '-----', &
                 partRow(4), '-----', partRow(5), '-----', partRow(6)

      if (i > 0) then
        write(*,'(A)') '|     |     |     |     |     |'
      end if
    end do
    write(*,*)
  end if

end subroutine showSimpleMeshPartitions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine visualizePartition()
  use mpi_h
  use zoltanRCBex
  implicit none

  !! Variables defined in zoltanRCBex module:
  !!      numLocObjs, GIDs, numExport, exportLocalGids, exportToPart

  !! local variables
  integer :: parts(numLocObjs)
  integer :: myrank, i, error

  call MPI_Comm_rank(MPI_COMM_WORLD, myrank, error)

  do i=1, numLocObjs, 1
    parts(i) = myRank;
  end do

  if (myRank== 0) then
    write (*,*) 'Mesh part assignments before calling Zoltan'
  end if

  call showSimpleMeshPartitions(myRank, numLocObjs, GIDs, parts);

  do i=1, numExport, 1
    parts(exportLocalGids(i)) = exportToPart(i)
  end do

  if (myRank == 0) then
    write (*,*) 'Mesh part assignments after calling Zoltan'
  end if

  call showSimpleMeshPartitions(myRank, numLocObjs, GIDs, parts)

end subroutine visualizePartition
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
program ofmd_fft
use, intrinsic :: iso_fortran_env, only: output_unit
use types, only: dp
use constants, only: i_, K2au, density2gcm3, u2au, s2au, Ha2eV
use md, only: velocity_verlet, positions_random, &
                calc_min_distance, positions_fcc, positions_bcc
use ewald_sums, only: ewald_box
use random, only: randn
use utils, only: init_random, stop_error, assert, linspace, clock, loadtxt
use ofdft, only: read_pseudo
use ofdft_fft, only: free_energy_min, radial_potential_fourier, &
    reciprocal_space_vectors, real2fourier, fourier2real, logging_info
use interp3d, only: trilinear
use poisson3d_assembly, only: func2quad
implicit none

! All variables are in Hartree atomic units

! XLBOMD parameters:
integer, parameter :: K = 5
real(dp), parameter :: kappa = 1.82_dp, alpha = 0.018_dp
integer, parameter :: c0 = -6, c1 = 14, c2 = -8, c3 = -3, c4 = 4, c5 = -1

integer :: N
integer :: i, j, steps, u
real(dp) :: dt, L, t, rho, scf_eps
real(dp), allocatable :: V(:, :), X(:, :), f(:, :), m(:)
real(dp), allocatable :: R(:), Ven_rad(:), &
    G(:, :, :, :), G2(:, :, :)
real(dp), allocatable :: Ven0G(:, :, :)
complex(dp), allocatable :: VenG(:, :, :), neG(:, :, :)
real(dp), allocatable :: ne(:, :, :), R2(:)
real(dp) :: Temp, Ekin, Temp_current, t3, t4
real(dp) :: Ediff, Z, Am
real(dp) :: Een_correction
real(dp) :: Eee, Een, Ts, Exc, Etot, Enn
real(dp), allocatable :: fnn(:, :), q(:), fen(:, :)
integer :: dynamics, functional, Ng, Nspecies, start, Nmesh
integer :: cg_iter
real(dp) :: sigma
real(dp), dimension(:,:,:,:), allocatable :: ne_aux
character(len=10) :: atom_name
real(dp) :: mdt1, mdt2

logging_info = .false. ! Turn of the INFO warnings

call read_input("OFMD.input", Temp, rho, Nspecies, N, Am, atom_name, start, &
    dynamics, functional, Ng, scf_eps, steps, dt)
call read_pseudo("fem/" // trim(atom_name) // ".pseudo", R, Ven_rad, Z, Ediff)
allocate(V(3, N), f(3, N), m(N))
allocate(Ven0G(Ng, Ng, Ng), VenG(Ng, Ng, Ng), ne(Ng, Ng, Ng), neG(Ng, Ng, Ng))
allocate(G(Ng, Ng, Ng, 3), G2(Ng, Ng, Ng))
allocate(fnn(3, N), q(N), fen(3, N))
allocate(ne_aux(Ng, Ng, Ng, -K:1))

q = Z
m = Am * u2au ! Using Hydrogen mass in atomic mass units [u]
L = (sum(m) / rho)**(1._dp/3)
print *, "----------------------------------------------------------------"
print *, "Input Summary:"
print *, "N =", N
print *, "Ng =", Ng
print *, "MD steps =", steps
print *, "rho = ",  rho, "a.u. =", rho * density2gcm3, "g/cm^3"
print *, "T =", Temp, "a.u. =", Temp / K2au, "K ="
print *, "  =", Temp * Ha2eV, "eV"
print "(' dt =', f8.2, ' a.u. = ', es10.2, ' s = ', es10.2, ' ps')", dt, &
    dt/s2au, dt/s2au * 1e12_dp
print "(' SCF_eps =', es10.2, ' a.u. = ', es10.2, ' eV')", scf_eps, &
    scf_eps * Ha2eV
print *, "Initial position:", start
print *, "Atomic mass:", Am, "u"
print *, "Atomic name: ", trim(atom_name)
print *, "----------------------------------------------------------------"
print *
print *, "Calculated quantities:"
print *, "L =", L, "a.u."
print *

print *, "Converting aperiodic radial Ven to periodic cartesian Ven"
call radial_potential_fourier(R, Ven_rad, L, Z, Ven0G, Een_correction)
print *, "  Done."
print *, "Atomic charge Z =", Z
Nmesh = 10000
allocate(R2(Nmesh))
R2 = linspace(0._dp, L/2, Nmesh)

call reciprocal_space_vectors([L, L, L], G, G2)

select case(start)
    case (0)
        print *, "Initial position: pos.txt"
        call loadtxt("pos.txt", X)
        call assert(size(X, 1) == 3)
        call assert(size(X, 2) == N)
    case (3)
        print *, "Initial position: FCC"
        allocate(X(3, N))
        call positions_fcc(X, L)
    case (4)
        print *, "Initial position: BCC"
        allocate(X(3, N))
        call positions_bcc(X, L)
    case (5)
        print *, "Initial position: random"
        ! Make it deterministic for now
        !call init_random()
        sigma = 1
        allocate(X(3, N))
        call positions_random(X, L, 2**(1._dp/6)*sigma, 10)
    case default
        call stop_error("Invalid initial condition.")
end select

print *, "Positions:"
do i = 1, N
    print *, i, X(:, i)
end do
print *, "Distances:"
do i = 2, N
    print *, i, calc_min_distance(X(:, :i-1), L, X(:, i))
end do
print *

! Initialize velocities based on Maxwell-Boltzmann distribution
call randn(V)
V = V * sqrt(Temp / spread(m, 1, 3))

Ekin = calc_Ekin(V, m)
Temp_current = 2*Ekin/(3*N)

! The average temperature (i.e. if we average Temp_current for many runs) will
! be Temp. But we want to set it exactly to Temp, so we rescale the velocities.
V = V * sqrt(Temp / Temp_current)

print *, "MD start:"

t = 0
ne = N * Z / L**3

open(newunit=u, file="ofmd_results.txt", status="replace")
call write_results_header(u)
close(u)

t = 0
call cpu_time(t3)
do i = 1, steps
    mdt1 = clock()
    print *, "Starting MD iteration:", i

    if (i > 10) then
        ! ne_aux(:, :, :,  1) ... n_{i+1}
        ! ne_aux(:, :, :,  0) ... n_{i}
        ! ne_aux(:, :, :, -1) ... n_{i-1}
        ! ne_aux(:, :, :, -2) ... n_{i-2}
        ! ne_aux(:, :, :, -3) ... n_{i-3}
        ! ne_aux(:, :, :, -4) ... n_{i-4}
        ! ne_aux(:, :, :, -5) ... n_{i-5}
        if (i < 100) then
            ne_aux(:, :, :, 1) = 2*ne - ne_aux(:, :, :, -1) &
                + alpha*(c0*ne &
                + c1*ne_aux(:, :, :, -1) + c2*ne_aux(:, :, :, -2) &
                + c3*ne_aux(:, :, :, -3) + c4*ne_aux(:, :, :, -4) &
                + c5*ne_aux(:, :, :, -5))
        else
            ne_aux(:, :, :, 1) = 2*ne_aux(:, :, :, 0) - ne_aux(:, :, :, -1) &
                + kappa*(ne-ne_aux(:, :, :, 0)) + alpha*(c0*ne_aux(:, :, :, 0) &
                + c1*ne_aux(:, :, :, -1) + c2*ne_aux(:, :, :, -2) &
                + c3*ne_aux(:, :, :, -3) + c4*ne_aux(:, :, :, -4) &
                + c5*ne_aux(:, :, :, -5))
        end if
    else
        ne_aux(:, :, :, 1) = ne
    end if

    do j = 0, -K, -1
        ne_aux(:, :, :, j) = ne_aux(:, :, :, j+1)
    end do

    ne = abs(ne_aux(:, :, :, 1))

    if (i == 1) then
        call forces(X, f)
    else
        call velocity_verlet(dt, m, L, forces, f, V, X)
    end if
    Ekin = calc_Ekin(V, m)
    Temp_current = 2*Ekin/(3*N)

    print *, "Nuclear forces:"
    print *, fnn(:, 1)
    print *, fnn(:, 2)
    print *, fnn(:, 3)
    print *, fnn(:, 4)

    print *, "Electronic forces:"
    print *, fen(:, 1)
    print *, fen(:, 2)
    print *, fen(:, 3)
    print *, fen(:, 4)

    print *, "total forces:"
    print *, f(:, 1)
    print *, f(:, 2)
    print *, f(:, 3)
    print *, f(:, 4)

    mdt2 = clock()
    open(newunit=u, file="ofmd_results.txt", position="append", status="old")
    call write_results_line(u)
    close(u)

    call write_results_header(output_unit)
    call write_results_line(output_unit)

    t = t + dt
end do
call cpu_time(t4)
close(u)
print "('TIMINGS')"
print "('MD run:      ', f10.4, 's')", t4-t3
print "('MD step:     ', f10.4, 's')", (t4-t3)/steps

contains

    subroutine forces(X, f)
    real(dp), intent(in) :: X(:, :) ! positions
    real(dp), intent(out) :: f(:, :) ! forces
    real(dp) :: stress(6)
    real(dp) :: fac(Ng, Ng, Ng)
    integer :: i
    print *, "-------------------------------------"
    print *, "Calculating forces"
    ! Calculate nuclear forces
    print *, "Calculating nuclear forces"
    call ewald_box(L, X, q, Enn, fnn, stress)

    ! Calculate the electronic forces
    print *, "Calculating VenG"
    VenG = 0
    do i = 1, N
        ! Note: we use minus sign in the forward real -> fourier transform.
        ! Then the following holds:
        !
        !         F[f(x+b)] = F[f(x)]*e^{+i*G*b},
        !
        ! with plus sign in the exponential. Finally, we are expressing
        ! Ven0(x-X) using the above formula (with b=-X) in terms of Ven0(x),
        ! i.e.:
        !
        !         F[Ven0(x-X)] = F[Ven0(x)]*e^{-i*G*X},
        !
        ! with minus sign in the exponential.
        VenG = VenG - Ven0G * exp(i_ * &
            (G(:,:,:,1)*X(1,i) + G(:,:,:,2)*X(2,i) + G(:,:,:,3)*X(3,i)))
    end do
    call assert(abs(VenG(1, 1, 1)) < epsilon(1._dp)) ! The G=0 component

    ! Energy calculation
    print *, "Minimizing free energy"
    call free_energy_min(N*Z, N, L, G2, Temp, VenG, ne, scf_eps, &
            Eee, Een, Ts, Exc, Etot, cg_iter)

    ! Forces calculation
    print *, "ne -> neG"
    call real2fourier(ne, neG)

    Een = Een + Een_correction * real(neG(1, 1, 1), dp) * N
    Etot = Eee + Een + Ts + Exc

    fen = 0
    print *, "Calculating fen"
    do i = 1, N
        ! We have minus sign in the exponential, per the definition of VenG
        ! above (see the comment there).
        fac = L**3*Ven0G*aimag(neG*exp(-i_ * &
            (G(:,:,:,1)*X(1,i) + G(:,:,:,2)*X(2,i) + G(:,:,:,3)*X(3,i))))
        fen(1, i) = sum(G(:,:,:,1)*fac)
        fen(2, i) = sum(G(:,:,:,2)*fac)
        fen(3, i) = sum(G(:,:,:,3)*fac)
    end do

    f = fnn + fen
    print *, "Done calculating forces."
    print *, "-------------------------------------"
    end subroutine

    real(dp) pure function calc_Ekin(V, m) result(Ekin)
    real(dp), intent(in) :: V(:, :) ! velocities
    real(dp), intent(in) :: m(:) ! masses
    Ekin = sum(m*sum(V**2, dim=1))/2
    end function

    subroutine read_input(filename, T, density, Nspecies, N, Am, Aname, &
        start, dynamics, functional, Ng, scf_eps, steps, dt)
    ! Reads the input file, returns values in a.u.
    character(len=*), intent(in) :: filename
    real(dp), intent(out) :: T, density, dt, scf_eps, Am
    character(len=*), intent(out) :: Aname
    integer, intent(out) :: Nspecies, N, start, dynamics, functional, Ng, steps
    integer :: AN
    real(dp) :: AZ
    integer :: u
    real(dp) :: skip
    open(newunit=u, file=filename, status="old")
    read(u, *) T, density, Nspecies, N
    read(u, *) start
    read(u, *) dynamics, functional, Ng, skip, skip, scf_eps
    read(u, *) steps, dt
    read(u, *) Aname, AZ, Am, AN
    close(u)
    ! Convert to atomic units:
    T = T / Ha2eV
    density = density / density2gcm3
    scf_eps = scf_eps / Ha2eV

    if (AN /= N) call stop_error("Inconsistent number of atoms")
    end subroutine

    subroutine write_results_header(u)
    integer, intent(in) :: u
    write(u, '(12a17)') "Time [a.u.]", "Fe [eV]", "Unn [eV]", "K [eV]", &
        "F [eV]", "T [eV]", "Ts [eV]", "Een [eV]", "Eee [eV]", "Exc [eV]", &
        "CG iter", "MD step time [s]"
    end subroutine

    subroutine write_results_line(u)
    integer, intent(in) :: u
    write(u, '(10f17.6, i17, f17.6)') &
        t, Etot*Ha2eV/N, Enn*Ha2eV/N, Ekin*Ha2eV/N, &
        (Etot + Enn + Ekin) * Ha2eV / N, Temp_current * Ha2eV, Ts * Ha2eV / N, &
        Een * Ha2eV / N, Eee * Ha2eV / N, Exc * Ha2eV / N, cg_iter, &
        mdt2 - mdt1
    end subroutine

end program

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
	PROGRAM FIXLIST
C
	CHARACTER*128 JUNK
C
	OPEN(UNIT=66,NAME='FFTMNB.LIS',STATUS='OLD')
	OPEN(UNIT=67,NAME='FFTMNB.FOR',STATUS='NEW')
C 1001	FORMAT(128A1)
 1001	FORMAT(A)
  100	READ(66,1001,END=200) JUNK
	WRITE(67,1001) JUNK(8:87)
	TYPE 1001, JUNK(8:87)
	GO TO 100
 200	CONTINUE
	CLOSE(UNIT=66)
	CLOSE(UNIT=67)
	STOP
	END

ext: for


text:
      subroutine pwrlat(word,wordp,lfno)
      use tfstk
      use ffs
      use sad_main
      use ffs_pointer, only:compelc
      use tffitcode
      implicit none
      type (sad_comp), pointer ::cmp
      integer*4 i,ln,lt,ni,nf,noel,l,ip,lfno,ielm,lene
      integer*4, parameter :: lline=131
      character*(*) word,wordp
      character*(MAXPNAME) name
      character*(lline) cline
      character patt*80
      logical psname,abbrev,tmatch,exist,patfit
      integer*4 ls(2),lf(2)
      data cline/' '/,ip/1/
c
      if(abbrev(word,'PS_NAME','_')) then
        psname=.true.
      else
        psname=.false.
      endif
      call getwdl2(word,wordp)
      if(word.eq.' ') then
        patt='*'
        go to 11
      endif
      do 10 i=1,nlat-1
        call compelc(i,cmp)
        if( tmatch(pname(cmp%id),wordp) ) then
          patt=wordp
          call getwdl2(word,wordp)
          goto 11
        endif
   10 continue
      patt=' '
   11 continue
      ni=ielm(wordp,exist)
      if(exist) then
        call getwdl2(word,wordp)
        nf=ielm(wordp,exist)
        if(exist) then
          call getwdl2(word,wordp)
        else
          nf=nlat
        endif
      else
        ni=1
        nf=nlat
      endif
c     print *,ni,nf
      if(ni.gt.nf) then
        ls(1)=ni
        lf(1)=nlat
        ls(2)=1
        lf(2)=nf
      else
        ls(1)=ni
        lf(1)=nf
        ls(2)=nlat
        lf(2)=1
      endif
      noel=0
      do 21 l=1,2
        do 20 i=ls(l),lf(l)
          call compelc(i,cmp)
          if(i .eq. nlat) then
            patfit='***'.eq.patt
     $           .or. patt .eq. '$$$' .or. patt.eq.'*'
          elseif(i .eq. 1)then
            patfit=tmatch(pname(cmp%id),patt) .or. patt .eq. '^^^'
          else
            patfit=tmatch(pname(cmp%id),patt)
          endif
          if( patfit ) then
            noel=noel+1
            if(i.eq.nlat) then
              name='$$$'
            else
              if(psname) then
                name=pname(cmp%id)
              else
                call elname(i,name)
              endif
            endif
            ln=len_trim(name)
            lt=(ln/8+1)*8
            if(ip+lt+1 .ge. lline)then
              write(lfno,'('' '',a)') cline
              cline=' '
              ip=1
            endif
            cline(ip:ip+lt-1)=name(1:ln)
            ip=ip+lt
          endif
   20   continue
   21 continue
      if(cline.ne.' ') then
        write(lfno,'('' '',a)') cline
        cline=' '
        ip=1
      endif
      write(lfno,'(I6,2A)') noel,' elements match ',patt(1:lene(patt))
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
      SUBROUTINE MRED2M (NUF,N2,NUS)
C
C     THIS SUBROUTINE FORMS THE HK MATRIX FOR THE MRED2 MODULE.
C
      INTEGER         DRY,GBUF1,OTFILE,Z,TYPIN,TYPOUT,HKPG,PHI12I,HK,
     1                RPRTN,CPRTN
      DIMENSION       ITRLR1(7),MODNAM(2),ITRLR2(7),ISUB(4),RZ(1)
      COMMON /BLANK / IDUM1,DRY,IDUM4,GBUF1,IDUM2(17),OTFILE(6),
     1                ISCR(10),KORLEN,KORBGN,IDUM3(14),NMODES
      COMMON /ZZZZZZ/ Z(1)
      COMMON /PACKX / TYPIN,TYPOUT,IROW,NROW,INCR
      EQUIVALENCE     (HK,ISCR(2)),(IDENT,ISCR(8)),(HKPG,ISCR(3)),
     1                (PHI12I,ISCR(8)),(CPRTN,ISCR(9)),(RPRTN,ISCR(9))
      EQUIVALENCE     (RZ(1),Z(1))
      DATA    MODNAM/ 4HMRED,4H2M  /
C
C     FORM HK MATRIX
C
C        **  **   **             **
C        *    *   *      .        *
C        * HK * = * HKPG . PHI12I *
C        *    *   *      .        *
C        **  **   **             **
C
      IF (DRY .EQ. -2) RETURN
      IF (NUF .EQ.  0) GO TO 30
      ITRLR1(1) = HKPG
      CALL RDTRL (ITRLR1)
      ITRLR2(1) = PHI12I
      CALL RDTRL (ITRLR2)
      INCR    = 1
      TYPIN   = 1
      TYPOUT  = 1
      IROW    = 1
      NROW    = ITRLR1(3) + ITRLR2(3)
      ISUB(1) = ITRLR1(3)
      ISUB(2) = ITRLR2(3)
      DO 20 I = 1,NROW
      RZ(KORBGN+I-1) = 0.0
      IF (I .GT. ITRLR1(3)) RZ(KORBGN+I-1) = 1.0
   20 CONTINUE
      IFORM = 7
      CALL MAKMCB (ITRLR2,RPRTN,NROW,IFORM,TYPIN)
      CALL GOPEN (RPRTN,Z(GBUF1),1)
      CALL PACK (Z(KORBGN),RPRTN,ITRLR2)
      CALL CLOSE (RPRTN,1)
      CALL WRTTRL (ITRLR2)
      ITYPE = 2
      CALL GMMERG (HK,HKPG,0,PHI12I,0,RPRTN,0,ISUB,ITYPE,Z(KORBGN),
     1            KORLEN)
      RETURN
C
C     NO UF POINTS
C
C        **  **   **     **
C        *    *   *   .   *
C        * HK * = * 0 . I *
C        *    *   *   .   *
C        **  **   **     **
C
   30 TYPIN  = 1
      TYPOUT = 1
      IROW   = 1
      NROW   = NMODES
      INCR   = 1
      IFORM  = 8
      IF (KORBGN+NMODES .GE. KORLEN) GO TO 100
C
C     GENERATE IDENTITY MATRIX
C
      CALL MAKMCB (ITRLR2,IDENT,NMODES,IFORM,TYPIN)
      CALL GOPEN  (IDENT,Z(GBUF1),1)
      DO 70 I = 1,NMODES
      DO 60 J = 1,NMODES
      RZ(KORBGN+J-1) = 0.0
      IF (J .EQ. I) RZ(KORBGN+J-1) = 1.0
   60 CONTINUE
   70 CALL PACK (Z(KORBGN),IDENT,ITRLR2)
      CALL CLOSE (IDENT,1)
      CALL WRTTRL (ITRLR2)
C
C     GENERATE ROW PARTITIONING VECTOR
C
      NROW = NUS + NMODES
      IF (KORBGN+NROW .GE. KORLEN) GO TO 100
      J = NROW
      DO 90 I = 1,J
      RZ(KORBGN+I-1) = 0.0
      IF (I .GT. NUS) RZ(KORBGN+I-1) = 1.0
   90 CONTINUE
      IFORM = 7
      CALL MAKMCB (ITRLR2,RPRTN,NROW,IFORM,TYPIN)
      CALL GOPEN (RPRTN,Z(GBUF1),1)
      CALL PACK (Z(KORBGN),RPRTN,ITRLR2)
      CALL CLOSE (RPRTN,1)
      CALL WRTTRL (ITRLR2)
C
C     FORM HK MATRIX
C
      ISUB(1) = NUS
      ISUB(2) = NMODES
      ITYPE   = 2
      CALL GMMERG (HK,0,0,IDENT,0,RPRTN,0,ISUB,ITYPE,Z(KORBGN),KORLEN)
      RETURN
C
C     PROCESS SYSTEM ERRORS
C
  100 IMSG  =-8
      IFILE = 0
      CALL SOFCLS
      CALL MESAGE (IMSG,IFILE,MODNAM)
      RETURN
C
      END


            program demo_vexit
            use M_draw, only: vinit, vexit, voutput, circle, linewidth, color
            use M_draw,    only  : D_BLACK,   D_WHITE
            use M_draw,    only  : D_RED,     D_GREEN,    D_BLUE
            use M_draw,    only  : D_YELLOW,  D_MAGENTA,  D_CYAN
            implicit none
            integer :: i

            ! assuming you have the NetPBM package installed
            ! set up output to create a GIF file called one.gif
            call voutput('|ppmtogif >vexit.3m_draw.gif')

            call vinit('p6') ! start graphics
            ! default window is -1 <= x <= 1, -1 <= y <= 1
            ! default viewport is left bottom justified square
            ! so essentially you have a square drawing surface
            ! with the point <0,0> at the center of the square
            ! with X and Y ranging from -1 to 1.

            call color(D_RED)
            call linewidth(100)
            ! this loop draws outside the current window
            ! but by default clipping is on
            do i=1,30
               call circle((real(i)-1)*0.04,0.0,1.0/real(i))
            enddo

            call vexit() ! exited graphics so can start again

            ! start second graphics session with new output
            ! device and output file
            call voutput('|ppmtogif >vexit2.3m_draw.gif')
            call vinit('p6')

            do i=10,1,-1
               call color(i)
               call circle(0.0,0.0,1.0/real(i))
            enddo

            call vexit()

            end program demo_vexit


text:
      REAL FUNCTION genexp(av)

C**********************************************************************
C
C     REAL FUNCTION GENEXP( AV )
C
C                    GENerate EXPonential random deviate
C
C
C                              Function
C
C
C     Generates a single random deviate from an exponential
C     distribution with mean AV.
C
C
C                              Arguments
C
C
C     AV --> The mean of the exponential distribution from which
C            a random deviate is to be generated.
C                              REAL AV
C
C     GENEXP <-- The random deviate.
C                              REAL GENEXP
C
C
C                              Method
C
C
C     Renames SEXPO from TOMS as slightly modified by BWB to use RANF
C     instead of SUNIF.
C
C     For details see:
C
C               Ahrens, J.H. and Dieter, U.
C               Computer Methods for Sampling From the
C               Exponential and Normal Distributions.
C               Comm. ACM, 15,10 (Oct. 1972), 873 - 882.
C
C**********************************************************************
C     .. Scalar Arguments ..
      REAL av
C     ..
C     .. External Functions ..
      REAL sexpo
      EXTERNAL sexpo
C     ..
C     .. Executable Statements ..
      genexp = sexpo()*av
      RETURN

      END

max_stars_repo_name: ppernot/MC-ChemDB

max_issues_repo_name: ppernot/MC-ChemDB

max_forks_repo_name: ppernot/MC-ChemDB


text:
Name: Susan Cameron

Website: http://www.des.ucdavis.edu/students/cameron/

Office: Wickson

Occupation: PhD candidate in http://ecology.ucdavis.edu Ecology.

Likes: the weekly veggie basket from the student farm, maps on computers, margaritas, bluegrass, Will Ferrell, and my hometown, Austin, Texas

Known for:  overstuffing burritos, drinking beer very slowly, saying interesting way too much, my obese, grumpy cat Oreo and the occasional bike accident. I have been the victim of many rutabaga breakins


text:
module probdata_module

      use amrex_fort_module, only : rt => amrex_real

!     These determine the refinement criteria
      real(rt), save ::    denerr,  dengrad
      real(rt), save ::    velerr,  velgrad
      real(rt), save ::  presserr,pressgrad
      real(rt), save ::   temperr, tempgrad
      real(rt), save ::    raderr,  radgrad
      integer , save ::  max_denerr_lev   ,max_dengrad_lev
      integer , save ::  max_velerr_lev   ,max_velgrad_lev
      integer , save ::  max_presserr_lev, max_pressgrad_lev
      integer , save ::  max_temperr_lev,  max_tempgrad_lev
      integer , save ::  max_raderr_lev,   max_radgrad_lev

      real(rt), save ::  center(3)

!     Sod variables
      real(rt), save ::  p_l, u_l, rho_l, p_r, u_r, rho_r, rhoe_l, rhoe_r, frac
      integer , save ::  nsub

!     These help specify which specific problem
      integer, save ::  probtype,idir

end module probdata_module

max_stars_repo_name: Alpine-DAV/Nyx

max_issues_repo_name: Alpine-DAV/Nyx

max_forks_repo_name: Alpine-DAV/Nyx