program data_misfit_02
    ! sum of misfits and check iterative status

    use seismo_parameters
    implicit none

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    integer, parameter :: NARGS = 5
    INTEGER :: isrc,iter
    INTEGER :: NPROC_DATA
    INTEGER :: i,j
    INTEGER :: ier
    real(kind=CUSTOM_REAL) :: misfit_cur
    real(kind=CUSTOM_REAL) :: misfit_initial
    character(len=MAX_STRING_LEN) :: arg(NARGS)
    character(len=MAX_STRING_LEN) :: directory
    character(len=MAX_FILENAME_LEN) :: FILENAME

    ! parse command line arguments
    if (command_argument_count() /= NARGS) then
        print *, 'USAGE: .bin/data_misfit.exe ...'
        stop ' Please check command line arguments'
    endif
    do i = 1, NARGS
    call get_command_argument(i,arg(i), status=ier)
    enddo
    read(arg(1),*) iter
    read(arg(2),*) step_length
    read(arg(3),*) compute_adjoint
    read(arg(4),*) NPROC_DATA
    directory=arg(5)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! sum event_misfit
    call sum_misfit(directory,misfit_cur,NPROC_DATA)
    !!! add current result to search history
    write(filename, "(a,i2.2)") trim(directory)//'/data_misfit/data_misfit_hist_iter',iter
    OPEN (IOUT, FILE=trim(filename),status='unknown',POSITION='APPEND')
    write(IOUT,'(f15.5,e15.5)') step_length,misfit_cur
    close(IOUT)

    ! check search status
    if (.not.compute_adjoint) then
        call check_linesearch(directory,iter)

    else
        if(iter==1)  then
            !!! misfit hist for iteration
            write(filename,'(a)') trim(directory)//'/data_misfit/data_misfit_hist.dat'
            OPEN (UNIT=IOUT, FILE=trim(filename),status='unknown',POSITION='APPEND')
            write(IOUT,'(I5,e15.5)') iter-1,misfit_cur
            close(IOUT)

            misfit_initial=misfit_cur
            !!! save initial absolute misfit
            write(filename,'(a)') trim(directory)//'/data_misfit/data_misfit_initial.dat'
            OPEN (UNIT=IOUT, FILE=trim(filename),status='unknown',POSITION='APPEND')
            write(IOUT,'(I5,e15.5)') iter-1,misfit_initial
            close(IOUT)

            write(filename,'(a)') trim(directory)//'/data_misfit/data_misfit_hist_norm.dat'
            OPEN (UNIT=IOUT, FILE=trim(filename),status='unknown',POSITION='APPEND')
            write(IOUT,'(I5,e15.5)') iter-1,misfit_cur/misfit_initial
            close(IOUT)
        endif
        ! search status initilization
        is_cont=1
        is_done=0
        is_brak=0
        next_step_length=initial_step_length_02
        optimal_step_length=0.0
        call check_iteration(directory)
    endif ! compute_adjoint

    !! SAVE search status
    write(filename,'(a)') trim(directory)//'/data_misfit/search_status.dat'
    OPEN (IOUT, FILE=trim(filename))
    write(IOUT,'(I5)') is_cont
    write(IOUT,'(I5)') is_done
    write(IOUT,'(I5)') is_brak
    write(IOUT,'(f15.5)') next_step_length
    write(IOUT,'(f15.5)') optimal_step_length
    close(IOUT)

end program data_misfit_02
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine sum_misfit(directory,misfit_cur,NPROC_DATA)
    use seismo_parameters
    implicit none
    integer :: ier,isrc,ip,NPROC_DATA
    real(kind=CUSTOM_REAL) :: temp,misfit_cur
    character(len=MAX_FILENAME_LEN) :: filename
    character(len=MAX_STRING_LEN) :: directory

    misfit_cur = 0.0_CUSTOM_REAL

    ! source loop
    do isrc=0, nsrc-1
    do ip=0,NPROC_DATA-1
    ! open file
    write(filename,'(a,i6.6,a,i6.6,a)') trim(directory)//'/',isrc,'/proc',ip,'_misfit.dat'
    OPEN (IIN,FILE= filename,STATUS='OLD',action='read',iostat=ier)
    if(ier>0) then
        print*,'Error opening event misfit file: ',filename
        stop
    else
        read(IIN,*) temp
    end if
    close(IIN)
    if(DISPLAY_DETAILS) then
        print*,'read misfit file from --- ',trim(filename), ' misfit=', temp
    endif
    !! sum over source (chi-squared)
    misfit_cur=misfit_cur+temp
    print*,'Summation of total data misfit becomes (PWY):', misfit_cur
    enddo !! ip
    enddo !! source loop
    print*,'Summation of total data misfit is (PWY):', misfit_cur
end subroutine sum_misfit
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine check_iteration(directory)
    use seismo_parameters
    implicit none
    integer ::  j,iter
    integer :: ier,i,niter
    real(kind=CUSTOM_REAL) :: temp(2)
    real(kind=CUSTOM_REAL) :: misfit_hist(iter_end)
    character(len=MAX_FILENAME_LEN) :: filename
    character(len=MAX_STRING_LEN) :: directory

    misfit_hist=0.0_CUSTOM_REAL

    write(filename, "(a)") trim(directory)//'/data_misfit/data_misfit_hist.dat'
    OPEN (IIN,FILE= filename,STATUS='OLD',action='read',iostat=ier)

    j=0
    do i=1,iter_end+1
    read(IIN,*,iostat=ier) iter,misfit_hist(iter+1)
    if (ier/=0) exit
    j=j+1
    enddo
    close(IIN)
    niter = j
    print*,'misfit_hist for niter ',niter, ' :',misfit_hist(1:niter)

    ! absolute misfit small enough
    if(misfit_hist(niter)<=SMALL_VAL) then
        is_cont=0
        is_brak=1
        next_step_length=0.0
        print*, 'stop due to current misfit value is small enough :',misfit_hist(niter)
    else if(j>1) then
        ! relative to initial misfit
        if(misfit_hist(niter)/misfit_hist(1) <=misfit_ratio_initial) then
            is_cont=0
            is_brak=1
            next_step_length=0.0
            print*, 'stop due to current misfit value is less than ',misfit_ratio_initial*100,'%',&
                ' relative to initial misfit :', &
                misfit_hist(niter)/misfit_hist(1)*100,'%'

            ! misfit reduction relative to previous iteration
        else if( (misfit_hist(niter-1) - misfit_hist(niter))/misfit_hist(niter-1) <=misfit_ratio_previous) then
            is_cont=0
            is_brak=0
            next_step_length=0.0
            print*, 'stop due to misfit reduction is less than ',misfit_ratio_initial*100,'%', &
                ' relative to previous iteration :',&
                (misfit_hist(niter-1) - misfit_hist(niter))/misfit_hist(niter-1) * 100, '%'
        endif
    endif

end subroutine check_iteration
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine check_linesearch(directory,iter)
    use seismo_parameters
    implicit none
    integer :: iter,j
    integer :: ier,i,step
    real(kind=CUSTOM_REAL) :: optimal_misfit
    real(kind=CUSTOM_REAL) :: temp(2)
    real(kind=CUSTOM_REAL) :: step_hist(max_step),misfit_hist(max_step)
    character(len=MAX_FILENAME_LEN) :: filename
    character(len=MAX_STRING_LEN) :: directory

    integer :: iter2
    integer :: ier2
    real(kind=CUSTOM_REAL) :: misfit_hist2(max_step)
    character(len=MAX_FILENAME_LEN) :: filename2

    step_hist = 0.0_CUSTOM_REAL
    misfit_hist=0.0_CUSTOM_REAL

    write(filename, "(a,i2.2)") trim(directory)//'/data_misfit/data_misfit_hist_iter',iter
    OPEN (IIN,FILE= filename,STATUS='OLD',action='read',iostat=ier)
    j=0
    do i=1,max_step+1
    read(IIN,*,iostat=ier) temp
    if (ier/=0) exit
    step_hist(i)=temp(1)
    misfit_hist(i)=temp(2)
    j=j+1
    enddo
    close(IIN)

    misfit_hist2=0.0_CUSTOM_REAL

    write(filename2, "(a)") trim(directory)//'/data_misfit/data_misfit_initial.dat'
    OPEN (IIN,FILE= filename2,STATUS='OLD',action='read',iostat=ier2)
    read(IIN,*,iostat=ier2) iter2,misfit_hist2(1)
    close(IIN)

    step=j

    !if(DISPLAY_DETAILS) then
    print*,'misfit_hist for nstep ',step
    print*,'step_hist : ',step_hist(1:step)
    print*,'misfit_hist: ',misfit_hist(1:step)
    !endif

    ! determine next step search status
    !  if(.not. backtracking) then
    print*,'line search method -- constant step size'

    if(step_length>=initial_step_length_02) then
        ! current status -- forward search
        if(misfit_hist(step)<misfit_hist(step-1))  then    ! decrease misfit
            ! two criteria : max_step; misfit reduction rate
            if(step<=max_step) then
                !.and.(misfit_hist(step-1) - misfit_hist(step))/misfit_hist(step-1) >= misfit_ratio_previous) then
                ! next status -- forward continue
                is_cont=1
                is_done=0
                is_brak=0
                next_step_length=step_hist(step)+initial_step_length_02
                optimal_step_length=step_hist(step)
                write(*,'(a,f15.5)') 'next step : forward continue -- next step_length=',next_step_length

            else
                is_cont=0
                is_done=1
                is_brak=0
                next_step_length=0.0
                optimal_step_length=step_hist(step)
                optimal_misfit=misfit_hist(step)
                write(*,'(a,f15.5)') 'next step : forward stop -- exceed max step, optimal step_length=',optimal_step_length
            endif

        else   ! not decrease misfit
            if(step_length>=1.5*initial_step_length_02) then
                !! more than one step forwad,  next status -- forward stop
                is_cont=0
                is_done=1
                is_brak=0
                next_step_length=0.0
                optimal_step_length=step_hist(step-1)
                optimal_misfit=misfit_hist(step-1)
                write(*,'(a,f15.5)') 'next step : forward done -- optimal step_length=',optimal_step_length

            else    !! next status -- backward start
                is_cont=1
                is_done=0
                is_brak=0
                next_step_length=step_hist(step)/2
                optimal_step_length=0.0
                write(*,'(a,f15.5)') 'next step : backward start -- next step_length=',next_step_length
            endif
        endif

    else   ! current status -- backward search
        if(misfit_hist(step)<misfit_hist(1)) then  ! next status -- backward stop
            is_cont=0
            is_done=1
            is_brak=0
            next_step_length=0.0
            optimal_step_length=step_hist(step)
            optimal_misfit=misfit_hist(step)
            write(*,'(a,f15.5)') 'next step : backward done -- optimal step_length=',optimal_step_length
        else
            if(step_length>min_step_length .and. step<=max_step) then !!  next status -- backward continue
                is_cont=1
                is_done=0
                is_brak=0
                next_step_length=step_hist(step)/2
                optimal_step_length=0.0
                write(*,'(a,f15.5)') 'next step : backward continue -- next step_length=',next_step_length

            else    !!  next status -- break
                is_cont=0
                is_done=0
                is_brak=1
                next_step_length=0.0
                optimal_step_length=0.0
                write(*,'(a)') 'next step : backward exit'
            endif
        endif
    endif

    if(is_done==1) then
        !!! misfit hist for iteration
        write(filename,'(a)') trim(directory)//'/data_misfit/data_misfit_hist.dat'
        OPEN (IOUT, FILE=filename,status='unknown',POSITION='APPEND')
        write(IOUT,'(I5,e15.5)') iter,optimal_misfit
        close(IOUT)

        write(filename,'(a)') trim(directory)//'/data_misfit/data_misfit_hist_norm.dat'
        OPEN (IOUT, FILE=filename,status='unknown',POSITION='APPEND')
        write(IOUT,'(I5,e15.5)') iter,optimal_misfit/misfit_hist2(1)
        close(IOUT)
        ! check iteration for next step
        call check_iteration(directory)
    endif
end subroutine check_linesearch


text:
module vd2d_mms_mod_2
    implicit none
    real, parameter :: pi = 3.1415926535897932
    real, private :: rho0,rho1,D,mu,uf,vf,a,k,b0,w,r10

    !functions
    PUBLIC :: vd2d_mms_init_2
    PUBLIC :: vd2d_mms_rho
    PUBLIC :: vd2d_mms_rho_of_z
    PUBLIC :: vd2d_mms_z
    PUBLIC :: vd2d_mms_z_src
    PUBLIC :: vd2d_mms_u
    PUBLIC :: vd2d_mms_u_src
    PUBLIC :: vd2d_mms_v
    PUBLIC :: vd2d_mms_v_src
    PUBLIC :: vd2d_mms_p
!————————————————————————————————————————————————————————————————————————————————————————
contains
!————————————————————————————————————————————————————————————————————————————————————————
    subroutine vd2d_mms_init_2 (den0,den1,diff,visc,xvel,yvel,ampl,wnum,delx,freq)
        implicit none
        real, intent (in) :: den0,den1,diff,visc
        real, intent (in) :: xvel,yvel,ampl,wnum,delx,freq
        rho0 = den0! density at z = 0
        rho1 = den1! density at z = 1
        mu = visc ! dynamic viscosity
        D = diff ! diffusion coefficient: rho * alpha_z
        uf = xvel ! x convective velocity
        vf = yvel ! y convective velocity
        a = ampl ! amplitude of rho perturbation
        k = wnum * pi! wave number of rho perturbation
        b0 = delx ! parameter controlling sharpness of transition
        w = freq ! inverse time-scale for ‘‘diffusion"
        r10 = rho1-rho0
    end subroutine vd2d_mms_init_2
    !————————————————————————————————————————————————————————————————————————————————————————
    function vd2d_mms_rho (x,y,t)
        implicit none
        real, intent (in) :: x,y,t
        real  :: vd2d_mms_rho
        real  :: z
        z = vd2d_mms_z (x,y,t)
        vd2d_mms_rho = vd2d_mms_rho_of_z (z)
    end function vd2d_mms_rho
    !————————————————————————————————————————————————————————————————————————————————————————
    function vd2d_mms_rho_of_z (z)
        implicit none
        real, intent (in) :: z
        real  :: vd2d_mms_rho_of_z
        vd2d_mms_rho_of_z = 1.0/(z/rho1+(1.0-z)/rho0)
    end function vd2d_mms_rho_of_z
    !————————————————————————————————————————————————————————————————————————————————————————
    function vd2d_mms_z (x,y,t)
        implicit none
        real, intent (in) :: x,y,t
        real  :: vd2d_mms_z
        real  :: xt,yt,xi,b,c
        xt = uf * t-x
        yt = vf * t-y
        b = b0 * exp (-w * t)
        xi = a * cos (k * yt)
        c = tanh (b * (xi + xt))
        vd2d_mms_z=(1.0 + c)/(1.0 + c+rho0/rho1 * (1.0-c))
    end function vd2d_mms_z
    !————————————————————————————————————————————————————————————————————————————————————————
    function vd2d_mms_z_src (x,y,t)
        implicit none
        real, intent (in) :: x,y,t
        real  :: vd2d_mms_z_src
        real  :: s0,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10, &
            s11,s12,s13,s14,s15,s16,s17,s18,s19,s20
        s0 = uf * t-x
        s1 = vf * t-y
        s2 = exp (-w * t)
        s3 = 1.0/s2
        s4 = cos (k * s1)
        s5 = sin (-k * s1)
        s6 = a * s4 + s0
        s7 = tanh (b0 * s2 * s6)
        s8 = exp (2.0 * b0 * s2 * s6)
        s9 = 1.0/2.0 + 1.0/2.0 * s7
        s10 = 1.0-s7 ** 2
        s11 = s8 + 1.0
        s12 = 2.0 * (w * a * s4 + w * uf * t-w * x-uf)
        s13 = s9 * r10 + rho0
        s14 = r10 * s7 + rho1 + rho0
        s15 = -s10 * b0 * s2 * r10
        s16 = -s10 * b0 * s2 * a * s5 * k * r10
        s17 = -b0 * w * s2 * s6 + b0 * s2 * (a * s5 * k * vf + uf)
        s18 = (1.0 + s7) * rho1
        s19 = 2.0 * (a * w * s4 * b0 + w * uf * t * b0-x * w * b0)-1.0/s11 * s12 * b0-w * s3 * log (s11)
        s20 = (-s11 ** 2 + s11 + s3 * s2 * s11 * s8) * w-s12 * b0 * s2 * s8
        vd2d_mms_z_src = &
            1.0/2.0 * s10 * s17 * r10 * s18/s14 + s13 * s10 * s17 * rho1/s14-s13 * s18/s14 ** 2 &
            * s10 * b0 * s2 * (-w * s6 + a * s5 * k * vf + uf) * r10 + 1.0/2.0 * s10 * s2 * rho1/s14 * r10 &
            * s19 + 1.0/2.0 * s18/s14 ** 2 * r10 * s19/b0 * s15-s18/s14 * r10 * s20/s11 ** 2 &
            -1.0/2.0 * s10 * b0 * s2 * a * s5 * k * r10 * s18/s14 * vf-s13 * s10 * b0 * s2 * a * s5 * k &
            * rho1/s14 * vf-s13 * s18/s14 ** 2 * vf * s16-D * (-2.0 * s7 * s10 * b0 ** 2 * s2 ** 2 * rho1 &
            /s14 + 2.0 * s10 * b0 * s2 * rho1/s14 ** 2 * s15 + 2.0 * s18/s14 ** 3 * s15 ** 2 + 2.0 &
            * s18/s14 ** 2 * s7 * s10 * b0 ** 2 * s2 ** 2 * r10-2.0 * s7 * s10 * b0 ** 2 * s2 ** 2 * a ** 2 * s5 &
            ** 2 * k ** 2 * rho1/s14-s10 * b0 * s2 * a * s4 * k ** 2 * rho1/s14 + 2.0 * s10 * b0 * s2 * a * s5 &
            * k * rho1/s14 ** 2 * s16 + 2.0 * s18/s14 ** 3 * s16 ** 2 + s18/s14 ** 2 * s10 * b0 * s2 * a * k &
            ** 2 * r10 * (2.0 * a * s7 * s5 ** 2 * s2 * b0 + s4))
    end function vd2d_mms_z_src
    !————————————————————————————————————————————————————————————————————————————————————————
    function vd2d_mms_u (x,y,t)
        implicit none
        real, intent (in) :: x,y,t
        real  :: vd2d_mms_u
        real  :: xt,yt,xi,b
        xt = uf * t-x
        yt = vf * t-y
        b = b0 * exp (-w * t)
        xi = a * cos (k * yt)
        vd2d_mms_u=&
            (rho1-rho0) * (-w * (xi + xt)+(w * (xi + xt)-uf)/(exp (2.0 * b * (xi + xt))+1.0) &
            + 1.0/2.0 * w/b * log (exp (2.0 * b * (xi + xt))+1.0)) &
            /vd2d_mms_rho (x,y,t)
    end function vd2d_mms_u
    !————————————————————————————————————————————————————————————————————————————————————————
    function vd2d_mms_u_src (x,y,t)
        implicit none
        real, intent (in) :: x,y,t
        real  :: vd2d_mms_u_src
        real  :: s0,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10, &
            s11,s12,s13,s14,s15,s16,s17,s18,s19,s20
        s0 = uf * t-x
        s1 = vf * t-y
        s2 = exp (-w * t)
        s3 = 1.0/s2
        s4 = cos (k * s1)
        s5 = sin (-k * s1)
        s6 = a * s4 + s0
        s7 = tanh (b0 * s2 * s6)
        s8 = exp (2.0 * b0 * s2 * s6)
        s9 = 1.0/2.0 + 1.0/2.0 * s7
        s10 = 1.0-s7 ** 2
        s11 = s8 + 1.0
        s12 = 2.0 * (w * a * s4 + w * uf * t-w * x-uf)
        s13 = s9 * r10 + rho0
        s14 = a * s5 * k * vf + uf
        s15 = log (s11) * s11
        s16 = r10 * (-w * a * s4-w * uf * t + w * x)+1.0/2.0 * r10/s11 * (s12 + w * s3 * s15/b0)
        s17=-k * a * s5 * r10 * (w * s11-s12 * b0 * s2 * s8 + s3 * w * s2 * s8 * s11-w * s11 ** 2)/s11 ** 2
        s18 = r10 * (-w * s11 ** 2-s12 * b0 * s2 * s8 + w * s11 + s3 * w * s2 * s8 * s11)/s11 ** 2
        s19 = s15 * w + 2.0 * (-b0 * s2 * s8 * s6 * w + b0 * s2 * s8 * a * s5 * k * vf + b0 * s2 * s8 * uf)
        s20=(s12 * b0 * s2 * s4 * s8 * s11 + w * s4 * s11 ** 3-w * s4 * s11 ** 2-2.0 * s12 * b0 ** 2 * s2 ** 2 &
            * a * s5 ** 2 * s8 * s11 + 2.0 * w * s3 * s2 ** 2 * a * s5 ** 2 * b0 * s8 * s11 ** 2-2.0 * w * s3 * s2 &
            ** 2 * a * s5 ** 2 * s8 ** 2 * b0 * s11 + 4.0 * s12 * b0 ** 2 * s2 ** 2 * a * s5 ** 2 * s8 ** 2-w * s2 * s3 &
            * s4 * s8 * s11 ** 2-4.0 * a * w * s5 ** 2 * b0 * s2 * s8 * s11) * a * k ** 2 * r10/s11 ** 3
        vd2d_mms_u_src=&
            -w * r10 * s14-1.0/2.0/s11 ** 2 * r10 * s12 * (-2.0 * b0 * w * s2 * s6 + 2.0 * b0 &
            * s2 * s14) * s8 + 1.0/2.0/s11 * r10 * (2.0 * w * a * s5 * k * vf + 2.0 * w * uf) &
            + 1.0/2.0 * s3 * w * r10 * s19/s11/b0-2.0 * s16/s13 * s18 + 1.0/2.0 &
            * s16 ** 2/s13 ** 2 * s10 * b0 * s2 * r10 + s17 * vf-mu * ((16.0/3.0/s11 ** 3 * r10 &
            * s12 * b0 ** 2 * s2 ** 2 * s8 ** 2-16.0/3.0/s11 ** 2 * r10 * w * b0 * s2 * s8-8.0 &
            /3.0/s11 ** 2 * r10 * s12 * b0 ** 2 * s2 ** 2 * s8-8.0/3.0 * s3 * w * s2 ** 2 * b0 * s8 &
            * r10 * (-s11 + s8)/s11 ** 2)/s13-4.0/3.0 * s18/s13 ** 2 * s10 * b0 * s2 * r10 &
            + 2.0/3.0 * s16/s13 ** 3 * s10 ** 2 * b0 ** 2 * s2 ** 2 * r10 ** 2 + 4.0/3.0 * s16 &
            /s13 ** 2 * s7 * s10 * b0 ** 2 * s2 ** 2 * r10 + s20/s13 + s17/s13 ** 2 * s10 * b0 * s2 * a * s5 * k &
            * r10 + 1.0/2.0 * s16/s13 ** 3 * s10 ** 2 * b0 ** 2 * s2 ** 2 * a ** 2 * s5 ** 2 * k ** 2 * r10 &
            ** 2 + s16/s13 ** 2 * s7 * s10 * b0 ** 2 * s2 ** 2 * a ** 2 * s5 ** 2 * k ** 2 * r10 + 1.0/2.0 &
            * s16/s13 ** 2 * s10 * b0 * s2 * a * s4 * k ** 2 * r10)
    end function vd2d_mms_u_src
    !————————————————————————————————————————————————————————————————————————————————————————
    function vd2d_mms_v (x,y,t)
        implicit none
        real, intent (in) :: x,y,t
        real  :: vd2d_mms_v
        !useless
        vd2d_mms_v = vf+x-x+y-y+t-t
    end function vd2d_mms_v
    !————————————————————————————————————————————————————————————————————————————————————————
    function vd2d_mms_v_src (x,y,t)
        implicit none
        real, intent (in) :: x,y,t
        real  :: vd2d_mms_v_src
        real  :: s0,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10, &
            s11,s12,s13,s14,s15,s16,s17,s18,s19
        s0 = uf * t-x
        s1 = vf * t-y
        s2 = exp (-w * t)
        s3 = 1.0/s2
        s4 = cos (k * s1)
        s5 = sin (-k * s1)
        s6 = a * s4 + s0
        s7 = tanh (b0 * s2 * s6)
        s8 = exp (2.0 * b0 * s2 * s6)
        s9 = 1.0/2.0 + 1.0/2.0 * s7
        s10 = 1.0-s7 ** 2
        s11 = s8 + 1.0
        s12 = 2.0 * (w * a * s4 + w * uf * t-w * x-uf)
        s13 = s9 * r10 + rho0
        s14=-r10 * (-w * s11 ** 2-s12 * b0 * s2 * s8 + w * s11 + s3 * w * s2 * s11 * s8)/s11 ** 2
        s15 = 2.0 * (-s12 * b0 * s2 * s8 + s11 * w)+s2 * s11 * (s12 * b0-s3 * w * s11 + s3 * w * s8)
        s16 = 2.0 * b0 * s11 * (w * a * s4 + w * uf * t-w * x)-s12 * b0-s3 * w * log (s11) * s11
        s17=-w * s11 ** 2-s12 * b0 * s2 * s8 + s11 * w + s3 * w * s2 * s11 * s8
        s18=-b0 * w * s2 * s6 + b0 * s2 * (a * s5 * k * vf + uf)
        s19 = 1.0/2.0 * (s10 * s18 * r10 * vf-s10 * b0 * s2 * a * s5 * k * r10 * vf ** 2)+s14 * vf
        vd2d_mms_v_src=&
            1.0/12.0 * mu * r10 ** 3 * s16 * b0/s11/s13 ** 3 * s10 ** 2 * s2 ** 2 * a * s5 * k &
            -mu * (1.0/6.0 * s14/s13 ** 2 * b0 * s2 * a * s5 * k * r10-1.0/6.0 * a &
            * s5 * k * r10 ** 2 * s17/s11 ** 2/s13 ** 2 * b0 * s2-1.0/6.0 * r10 ** 2 * s16 &
            * b0/s11/s13 ** 2 * s7 * s2 ** 2 * a * s5 * k) * s10 + s19 + 2.0/3.0 * mu * b0 * s2 &
            * s8 * a * s5 * k * r10 * s15/s11 ** 3/s13
    end function vd2d_mms_v_src
    !————————————————————————————————————————————————————————————————————————————————————————
    function vd2d_mms_p (x,y,t)
        implicit none
        real, intent (in) :: x,y,t
        real  :: vd2d_mms_p
        vd2d_mms_p = 0.0*x*y*t
    !! NOTE: pressure can vary by a constant and still satisfy the mms
    end function vd2d_mms_p

!useless functions here
!subroutine vd1d_mms_init (rho_0,rho_1,k_1,k_2,width,diff)
!implicit none
!real, intent (in) :: rho_0,rho_1,k_1,k_2,width,diff
!k1 = k_1 ! inverse time-scale for amplitude decay
!k2 = k_2 ! inverse time-scale for width decay
!rho0 = rho_0! density at z = 0
!rho1 = rho_1! density at z = 1
!w0 = width ! initial width
!D = diff ! z-diffusivity
!r01 = rho0-rho1
!k12 = k1-k2
!k21 = k2-k1   !diff
!
end module vd2d_mms_mod_2

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


text:
      SUBROUTINE MULP(XXXX,YYYY,ZZZZ)
C
C THIS SUBROUTINE PERFORMS MULTIPLICATION OF COMPLEX NUMBERS:
C ZZZZ = XXXX*YYYY
C
C NOTE:  IN THE CALLING ROUTINE, ZZZZ SHOULD NOT BE THE SAME
C AS XXXX OR YYYY.  HOWEVER, XXXX MAY BE THE SAME AS YYYY.
C THUS, "CALL MULP(X,X,Z)" IS OK, BUT "CALL MULP(X,Y,X)" IS NOT.
C
C ON INPUT:
C
C XXXX  IS AN ARRAY OF LENGTH TWO REPRESENTING THE FIRST COMPLEX
C       NUMBER, WHERE XXXX(1) = REAL PART OF XXXX AND XXXX(2) =
C       IMAGINARY PART OF XXXX.
C
C YYYY  IS AN ARRAY OF LENGTH TWO REPRESENTING THE SECOND COMPLEX
C       NUMBER, WHERE YYYY(1) = REAL PART OF YYYY AND YYYY(2) =
C       IMAGINARY PART OF YYYY.
C
C ON OUTPUT:
C
C ZZZZ  IS AN ARRAY OF LENGTH TWO REPRESENTING THE RESULT OF
C       THE MULTIPLICATION, ZZZZ = XXXX*YYYY, WHERE ZZZZ(1) =
C       REAL PART OF ZZZZ AND ZZZZ(2) = IMAGINARY PART OF ZZZZ.
C
C DECLARATION OF INPUT
      DOUBLE PRECISION XXXX,YYYY
      DIMENSION XXXX(2),YYYY(2)
C
C DECLARATION OF OUTPUT
      DOUBLE PRECISION ZZZZ
      DIMENSION ZZZZ(2)
C
      ZZZZ(1) = XXXX(1)*YYYY(1) - XXXX(2)*YYYY(2)
      ZZZZ(2) = XXXX(1)*YYYY(2) + XXXX(2)*YYYY(1)
      RETURN
      END

ext: f

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
submodule (io:plasma_input) plasma_input_raw

use timeutils, only : date_filename

implicit none (type, external)

contains

module procedure input_root_currents_raw
  !! READS, AS INPUT, A FILE GENERATED BY THE GEMINI.F90 PROGRAM

  real(wp), dimension(:,:,:), allocatable :: tmparray3D
  real(wp), dimension(:,:,:,:), allocatable :: tmparray4D
  character(:), allocatable :: filenamefull
  real(wp), dimension(:,:,:), allocatable :: J1all,J2all,J3all
  real(wp), dimension(:,:,:), allocatable :: tmpswap
  real(wp) :: tmpdate

  !>  CHECK TO MAKE SURE WE ACTUALLY HAVE THE DATA WE NEED TO DO THE MAG COMPUTATIONS.
  if (flagoutput==3) error stop '  !!!I need current densities in the output to compute magnetic fields!!!'

  !> FORM THE INPUT FILE NAME
  filenamefull = date_filename(outdir,ymd,UTsec) // '.dat'
  print *, 'Input file name for current densities:  ', filenamefull

  block
    integer :: u
    open(newunit=u,file=filenamefull,status='old',form='unformatted',access='stream',action='read')
    read(u) tmpdate
    print *, 'File year:  ',tmpdate
    read(u) tmpdate
    print *, 'File month:  ',tmpdate
    read(u) tmpdate
    print *, 'File day:  ',tmpdate
    read(u) tmpdate
    print *, 'File UThrs:  ',tmpdate

    !> LOAD THE DATA
    if (flagoutput==2) then    !the simulation data have only averaged plasma parameters
      print *, '  Reading in files containing averaged plasma parameters of size:  ',lx1*lx2all*lx3all
      allocate(tmparray3D(lx1,lx2all,lx3all))
      !MZ:  I've found what I'd consider to be a gfortran bug here.  If I read
      !in a flat array (i.e. a 1D set of data) I hit EOF, according to runtime
      !error, well before I'm actually out of data this happens with a 20GB
      !input file for not for a 3GB input file...  This doesn't happen when I do
      !the reading with 3D arrays.
      read(u) tmparray3D    !ne - could be done with some judicious fseeking...
      read(u) tmparray3D    !vi
      read(u) tmparray3D    !Ti
      read(u) tmparray3D    !Te
      deallocate(tmparray3D)
    else    !full output parameters are in the output files
      print *, '  Reading in files containing full plasma parameters of size:  ',lx1*lx2all*lx3all*lsp
      allocate(tmparray4D(lx1,lx2all,lx3all,lsp))
      read(u) tmparray4D
      read(u) tmparray4D
      read(u) tmparray4D
      deallocate(tmparray4D)
    end if

    !> PERMUTE THE ARRAYS IF NECESSARY
    print *, '  File fast-forward done, now reading currents...'
    allocate(J1all(-1:lx1+2,-1:lx2all+2,-1:lx3all+2),J2all(-1:lx1+2,-1:lx2all+2,-1:lx3all+2), &
               J3all(-1:lx1+2,-1:lx2all+2,-1:lx3all+2))
    !! no need to permute dimensions for 3D simulations
    read(u) J1all(1:lx1,1:lx2all,1:lx3all),J2all(1:lx1,1:lx2all,1:lx3all),J3all(1:lx1,1:lx2all,1:lx3all)
    close(u)
  end block
  print *, 'Min/max current data:  ',minval(J1all(1:lx1,1:lx2all,1:lx3all)),maxval(J1all(1:lx1,1:lx2all,1:lx3all)), &
                                     minval(J2all(1:lx1,1:lx2all,1:lx3all)),maxval(J2all(1:lx1,1:lx2all,1:lx3all)), &
                                     minval(J3all(1:lx1,1:lx2all,1:lx3all)),maxval(J3all(1:lx1,1:lx2all,1:lx3all))

  !> DISTRIBUTE DATA TO WORKERS AND TAKE A PIECE FOR ROOT
  call bcast_send3D_ghost(J1all,tag%J1,J1)
  call bcast_send3D_ghost(J2all,tag%J2,J2)
  call bcast_send3D_ghost(J3all,tag%J3,J3)

  !> CLEAN UP MEMORY
  deallocate(J1all,J2all,J3all)
end procedure input_root_currents_raw

module procedure input_root_mpi_raw
  !! READ INPUT FROM FILE AND DISTRIBUTE TO WORKERS.
  !! STATE VARS ARE EXPECTED INCLUDE GHOST CELLS.  NOTE ALSO
  !! THAT RECORD-BASED INPUT IS USED SO NO FILES > 2GB DUE
  !! TO GFORTRAN BUG WHICH DISALLOWS 8 BYTE INTEGER RECORD
  !! LENGTHS.

  integer :: lx1,lx2,lx3,lx2all,lx3all,isp

  real(wp), dimension(-1:size(x1,1)-2,-1:size(x2all,1)-2,-1:size(x3all,1)-2,1:lsp) :: nsall, vs1all, Tsall
  real(wp), dimension(:,:,:,:), allocatable :: statetmp
  integer :: lx1in,lx2in,lx3in,u, utrace
  real(wp) :: tin

  real(wp) :: tstart,tfin

  !> so that random values (including NaN) don't show up in Ghost cells
  nsall = 0._wp
  vs1all= 0._wp
  Tsall = 0._wp
  ns=0._wp
  vs1=0._wp
  Ts=0._wp

  !> SYSTEM SIZES
  lx1=size(ns,1)-4
  lx2=size(ns,2)-4
  lx3=size(ns,3)-4
  lx2all=size(x2all)-4
  lx3all=size(x3all)-4

  !> READ IN FROM FILE, AS OF CURVILINEAR BRANCH THIS IS NOW THE ONLY INPUT OPTION
  call get_simsize3(indatsize, lx1in, lx2in, lx3in)
  print *, 'Input file has size:  ',lx1in,lx2in,lx3in
  print *, 'Target grid structure has size',lx1,lx2all,lx3all,lx2,lx3

  if (.not. (lx1==lx1in .and. lx2all==lx2in .and. lx3all==lx3in)) then
    error stop 'The input data must be the same size as the grid which you are running the simulation on' // &
         '- use a script to interpolate up/down to the simulation grid'
  end if

  block
  integer :: u
  real(wp), dimension(3) :: ymdtmp
  open(newunit=u,file=indatfile,status='old',form='unformatted', access='stream', action='read')
  read(u) ymdtmp,tin

  read(u) nsall(1:lx1,1:lx2all,1:lx3all,1:lsp)
  read(u) vs1all(1:lx1,1:lx2all,1:lx3all,1:lsp)
  read(u) Tsall(1:lx1,1:lx2all,1:lx3all,1:lsp)
  close(u)
  end block

  !> Don't support restarting potential right now for raw input
  Phiall=0._wp
  Phi=0._wp

  !> Check the parameters (full grid)
  !  print*, 'Full grid parameter synopsis:  '
  !  print '(A,2ES11.2)', 'Min/max input density:',     minval(nsall(:,:,:,7)),  maxval(nsall(:,:,:,7))
  !  print '(A,2ES11.2)', 'Min/max input velocity:',    minval(vs1all(:,:,:,:)), maxval(vs1all(:,:,:,:))
  !  print '(A,2ES11.2)', 'Min/max input temperature:', minval(Tsall(:,:,:,:)),  maxval(Tsall(:,:,:,:))

  !> ROOT BROADCASTS IC DATA TO WORKERS
  call cpu_time(tstart)
  call bcast_send(nsall,tag%ns,ns)
  call bcast_send(vs1all,tag%vs1,vs1)
  call bcast_send(Tsall,tag%Ts,Ts)
  call bcast_send3D_ghost(Phiall,tag%Phi,Phi)
  call cpu_time(tfin)
  print *, 'Done sending ICs to workers...  CPU elapsed time:  ',tfin-tstart

  !> Check the parameters (root)
  !  print*, 'Root chunk parameter synopsis:  '
  !  print '(A,2ES11.2)', 'Min/max input density:',     minval(ns(:,:,:,7)),  maxval(ns(:,:,:,7))
  !  print '(A,2ES11.2)', 'Min/max input velocity:',    minval(vs1(:,:,:,:)), maxval(vs1(:,:,:,:))
  !  print '(A,2ES11.2)', 'Min/max input temperature:', minval(Ts(:,:,:,:)),  maxval(Ts(:,:,:,:))
end procedure input_root_mpi_raw

end submodule plasma_input_raw

max_forks_repo_name: mattzett/GEMINI


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

C=======================================================================
      LOGICAL FUNCTION CKCROS (IMID, IH, NLNKF, LINKF1, XN, YN, ZN)
C=======================================================================

C   --*** CKCROS *** (MESH) Hide partial lines by face
C   --   Written by Amy Gilkey - revised 10/02/86
C   --
C   --CKCROS deletes partial lines that have a visible node along the edge
C   --of a quarilateral and the hidden node is within the area defined
C   --by the vectors from the corner to the two adjacent nodes.  Such
C   --a face totally hides the partial line.
C   --
C   --Parameters:
C   --   IMID     - IN - the node in the face that is closest to the
C   --                   visible node
C   --   IH       - IN - the hidden nodes
C   --   NLNKF    - IN - the number of nodes per face
C   --   LINKF1   - IN - the connectivity for the face containing the edge
C   --   XN,YN,ZN - IN - the nodal coordinates

      include 'debug.blk'
C     common /debugc/ cdebug
C     common /debugn/ idebug
C     character*8 cdebug

      INTEGER LINKF1(NLNKF)
      REAL XN(*), YN(*), ZN(*)

      CKCROS = .FALSE.

C   --Find left and right node for face

      DO 100 ILINK = 1, 4
         IF (LINKF1(ILINK) .EQ. IMID) IXMID = ILINK
  100 CONTINUE
      IF (IXMID .LT. 4) THEN
         ILFT = LINKF1(IXMID+1)
      ELSE
         ILFT = LINKF1(1)
      END IF
      IF (IXMID .GT. 1) THEN
         IRGT = LINKF1(IXMID-1)
      ELSE
         IRGT = LINKF1(4)
      END IF
      IF ((ILFT .EQ. IH) .OR. (IRGT .EQ. IH)) GOTO 110
      XMID = XN(IMID)
      YMID = YN(IMID)
      XLFT = XN(ILFT)
      YLFT = YN(ILFT)
      XRGT = XN(IRGT)
      YRGT = YN(IRGT)
      IF ((XLFT-XMID)*YRGT + (XRGT-XLFT)*YMID
     &   + (XMID-XRGT)*YLFT .LT. 0.0) THEN
         I = IRGT
         IRGT = ILFT
         ILFT = I
         XLFT = XN(ILFT)
         YLFT = YN(ILFT)
         XRGT = XN(IRGT)
         YRGT = YN(IRGT)
      END IF

C   --Check 2 triangles formed by 2 adjacent nodes of the face
C   --and the node have a positive area

      X0 = XN(IH)
      Y0 = YN(IH)
      ALFT = (XLFT-XMID)*Y0 + (X0-XLFT)*YMID + (XMID-X0)*YLFT
      IF (ALFT .LT. 0.0) GOTO 110
      ARGT = (XMID-XRGT)*Y0 + (X0-XMID)*YRGT + (XRGT-X0)*YMID
      IF (ARGT .LT. 0.0) GOTO 110

C   --Check if the node's Z is behind the face's Z at that point

      ZMIN = MIN
     &   (ZN(LINKF1(1)), ZN(LINKF1(2)), ZN(LINKF1(3)), ZN(LINKF1(4)))

      IF (ZMIN .LE. ZN(IH)) THEN

C      --Calculate normal at visible node with diagonal node and
C      --left or right node

         IF (IXMID .GT. 2) THEN
            IDIA = LINKF1(IXMID-2)
         ELSE
            IDIA = LINKF1(IXMID+2)
         END IF

         ZMID = ZN(IMID)

         AMIN = MIN (ALFT, ARGT)

         IF (AMIN .EQ. ALFT) THEN
            ISIDA = ILFT
            ISIDB = IDIA
         ELSE IF (AMIN .EQ. ARGT) THEN
            ISIDA = IDIA
            ISIDB = IRGT
         END IF
         AX = XN(ISIDA) - XMID
         AY = YN(ISIDA) - YMID
         AZ = ZN(ISIDA) - ZMID
         BX = XN(ISIDB) - XMID
         BY = YN(ISIDB) - YMID
         BZ = ZN(ISIDB) - ZMID
         XNORM = 0.5 * (AY*BZ - BY*AZ)
         YNORM = 0.5 * (AZ*BX - BZ*AX)
         ZNORM = 0.5 * (AX*BY - BX*AY)

         IF ((XNORM*(X0-XMID) + YNORM*(Y0-YMID)
     &      + ZNORM*(ZN(IH)-ZMID)) .GE. 0.0) GOTO 110
      END IF

C   --Mark totally hidden lines

      CKCROS = .TRUE.

  110 CONTINUE
      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
        interface
           subroutine sub(char1,char2,vec)
     	   character(kind=1,len=7),dimension(:,:,:) :: char1,char2
     	   integer(kind=4),dimension(:) :: vec
           end subroutine
        end interface
	character(kind=1,len=7),dimension(6,6,6) :: ch1,ch2
	integer(kind=4),dimension(6)  :: v
	integer(kind=4)               :: i,j,k
        do i=1,6
           do j=1,6
              do k=1,6
                 ch1(i,j,k) = 1_' '
                 ch2(i,j,k) = 1_' '
              end do
           end do
           v(i) = i
        end do
	end program

	subroutine sub(ch1,ch2,v)
	character(kind=1,len=7),dimension(:,:,:) :: ch1,ch2
	integer(kind=4)::i,j,k
	integer(kind=4),dimension(:) :: v
        logical(kind=1)::flag
	end subroutine

max_stars_repo_name: xevolver/xevxml

max_issues_repo_name: xevolver/xevxml

max_forks_repo_name: xevolver/xevxml


text:
module cuda_pinned_builder_module
    use :: memory_api, only : data_storage
    use :: host_builder_module, only : host_builder
    use :: cuda_pinned_storage_module, only : cuda_pinned_storage

    implicit none
    private

    public :: cuda_pinned_builder

    type, extends(host_builder) :: cuda_pinned_builder
    contains
        procedure :: allocate_class => allocate_class
    end type cuda_pinned_builder

contains
    subroutine allocate_class(this, dst)
        class(cuda_pinned_builder), intent(in) :: this
        class(data_storage), allocatable, intent(out) :: dst

        allocate(cuda_pinned_storage::dst)
        call dst%clear()
    end subroutine allocate_class
end module cuda_pinned_builder_module

max_stars_repo_name: cheshyre/ntcl-data

max_issues_repo_name: cheshyre/ntcl-data

max_forks_repo_name: cheshyre/ntcl-data


text:
program clauses

   implicit none
   integer, parameter :: n=10
   integer :: i, a, b, c, d
   character(len=*), parameter :: layout_1 = "(a, 4(a, i12))"
   character(len=*), parameter :: layout_2 = "(a, i4, 4(a, i12))"

   a = 1
   b = 1
   c = 1
   d = 1

   print layout_1, "BEFORE", " SHARED:", a, " PRIVATE:", b, " FIRSTPRIVATE:", c, " LASTPRIVATE:", d

!$OMP PARALLEL DO SHARED(a) PRIVATE(b) FIRSTPRIVATE(c) LASTPRIVATE(d)
   do i = 1, n

      if (n<11) print layout_2, "I=", i, " SHARED:", a, " PRIVATE:", b, " FIRSTPRIVATE:", c, " LASTPRIVATE:", d
      a = a + i
      b = b + i
      c = c + i
      d = d + i
   end do
!$OMP END PARALLEL DO

   print layout_1, "AFTER ", " SHARED:", a, " PRIVATE:", b, " FIRSTPRIVATE:", c, " LASTPRIVATE:", d

end program clauses

max_stars_repo_name: WVUHPC/Modern-Fortran

max_issues_repo_name: WVUHPC/Modern-Fortran

max_forks_repo_name: WVUHPC/Modern-Fortran


text:
| demonstrate bug [3a39ec2726d36e09b414d5bad1e21abb7e61326e]

:   buf:  create  0 , 60 allot  does> ;

buf: tmp  | the following is ignored:

." You're lucky if you see this." cr bye

max_stars_repo_name: ronaaron/reva

max_issues_repo_name: ronaaron/reva

max_forks_repo_name: ronaaron/reva


text:
************************************************************************
*
*     Function LOGLK                                  Called by: CALIBR
*
*     calculate the log likelihood of Tobit regression (Paarsch, 1984)
*
*     locals vars
*     -----------
*     YHAT    predicted value of the response variable
*
************************************************************************
      DOUBLE PRECISION FUNCTION LOGLK(NPAR,NOBSC,XLCAL,YLCAL,CENSFLAG,
     &                                PARAML,RESSTD,NOBSCI)
*
*     subroutine arguments
*
      INTEGER*4 NPAR,NOBSC,NOBSCI
      LOGICAL CENSFLAG(*)
      DOUBLE PRECISION RESSTD,YLCAL(*),PARAML(*),XLCAL(NOBSCI,*)
*
*     local vars
*
      INTEGER*4 I
      DOUBLE PRECISION YHAT
*
*     function declarations
*
      DOUBLE PRECISION ERF,PRED
*
*     compute the log likelihood due to censored and uncensored data
*
*     if the observation is censored, its contribution to LOGLK is based
*     on the integral of the probability density of a standardized
*     normal variable from X to infinity (the complement of the
*     cumulative distribution function)(Abramowitz and Stegun, 1964, p.
*     297 & 931). If the observation is uncensored, the contribution is
*     based on the probability density function of the standard normal
*     distribution, divided by the residual standard deviation (RESSTD).
*
      LOGLK = 0.D0

      DO 10 I=1,NOBSC
         YHAT = PRED(NPAR,NOBSCI,I,XLCAL,PARAML)
         IF (CENSFLAG(I)) THEN
            LOGLK = LOGLK
     &             + DLOG((1.D0
     &                -ERF((-(YLCAL(I)-YHAT)/RESSTD)/DSQRT(2.D0)))/2.D0)
         ELSE
            LOGLK = LOGLK
     &             + DLOG(((1.D0/DSQRT(2.D0*3.14159265359D0)) *
     &                     DEXP(-(((YLCAL(I)-YHAT)/RESSTD)**2)/2.D0)
     &                                                     /RESSTD))
         ENDIF
 10   CONTINUE

      RETURN
      END

max_stars_repo_name: lthiamodelers/baseflow-coefficients

max_issues_repo_name: lthiamodelers/baseflow-coefficients

max_forks_repo_name: lthiamodelers/baseflow-coefficients


text:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Compute the matrix R.
c This version detrends every search window W before computing rho(A,W).

      subroutine compute_R_detrend(A,ia,ja,B,ib,jb,R,ir,jr,
     &		W,rsw,csw)
      implicit none
c Inputs...
      integer ia,ja
      integer*2 A(-ia:ia,-ja:ja)	!nlw x nsw shape function.
      integer ib,jb
      integer*2 B(-ib:ib,-jb:jb)	!search window centered at D(lnom,snom).
      integer ir,jr
c Scratch space...
      integer*2 W(-ia:ia,-ja:ja)
      integer rsw(-ja:ja)		!sums along each row of w
      integer csw(-ib:ib)		!sums along each column of B
c Output...
      real*4 R(-ir:ir,-jr:jr)		!rho for nver x nhor search area

c Template A mean and sigma...
      common/ca/area,meanA,sigA
      real area,meanA,sigA

c Local variables...
      integer ii,jj		!search window offsets
      integer i,j		!A(i,j) = B(i+ii,j+jj)

      integer sw,sw2,saw
      real meanW,varW,covAW,sigW

      do 20 jj=-jr,jr		!row loop thru R
      do 10 ii=-ir,ir		!column loop thru R
      call detrend_window(B,ib,jb,ii,jj,W,ia,ja,rsw,csw)
      sw = 0
      sw2 = 0
      saw = 0
      do j=-ja,ja
         do i=-ia,ia
            sw = sw + w(i,j)
            sw2 = sw2 + w(i,j)**2
            saw = saw + a(i,j)*w(i,j)
         enddo
      enddo
      meanW = sw/area
      varW = sw2/area - meanW**2
      if (varW.gt.0) then	!do this if there is a signal
         covAW = saw/area - meanA*meanW
         sigW = sqrt(varW)
         R(ii,jj) = covAW/(sigA*sigW)
      else
         R(ii,jj) = -999	!no pulse, rho has no meaning.
      endif
   10 continue
   20 continue

      return
      end

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
MODULE maskf_module

  implicit none

  PRIVATE
  PUBLIC :: nmsk,maskr,xm,dxm,makemaskf

  integer,parameter :: nmsk=201
  real(8) :: maskr(0:nmsk),xm(0:nmsk),dxm

CONTAINS

  SUBROUTINE makemaskf(eta)
    real(8),intent(IN) :: eta
    real(8),allocatable :: X12(:,:),iX12(:,:),lambda(:),rwork(:)
    real(8) :: dx,x1,x2,x,c,pi
    integer,parameter :: matz=1
    integer :: ierr,mxtmp,i,j,LWORK,LIWORK
    integer,allocatable :: iwork(:)
    pi=acos(-1.d0)
    mxtmp=nmsk-1
    dx=1.d0/mxtmp
    allocate( lambda(mxtmp) )
    allocate( X12(mxtmp,mxtmp) )
    do j=1,mxtmp
       x2=j*dx
       do i=j,mxtmp
          x1=i*dx
          if ( i==j ) then
             X12(i,j)=sin((x1+x2)*eta)/(x1+x2)+mxtmp*pi-eta
          else
             X12(i,j)=sin((x1+x2)*eta)/(x1+x2)-sin((x1-x2)*eta)/(x1-x2)
             X12(j,i)=X12(i,j)
          end if
       end do
    end do
    LWORK=1+6*mxtmp+2*mxtmp*mxtmp
    LIWORK=3+5*mxtmp
    allocate( rwork(LWORK),iwork(LIWORK) )
    call DSYEVD('V','L',mxtmp,X12,mxtmp,lambda,rwork,LWORK,iwork,LIWORK,ierr)
    deallocate( iwork,rwork )
    c=X12(1,1)/dx
    do i=1,mxtmp
       x=i*dx
       maskr(i)=X12(i,1)/x/c
    end do
    maskr(0)=1.d0
    maskr(nmsk)=0.d0
    dxm=dx
    do i=0,nmsk
       xm(i)=i*dxm
    end do
    deallocate( lambda )
    deallocate( X12 )
    return
  END SUBROUTINE makemaskf

END MODULE maskf_module

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
! <cphotr.f90 - A component of the City-scale
!                 Chemistry Transport Model EPISODE-CityChem>
!*****************************************************************************!
!*
!* EPISODE - An urban-scale air quality model
!* ==========================================
!* Copyright (C) 2018  NILU - Norwegian Institute for Air Research
!*                     Instituttveien 18
!*                     PO Box 100
!*                     NO-2027 Kjeller
!*                     Norway
!*
!*                     Contact persons: Gabriela Sousa Santos - gss@nilu.no
!*                                      Paul Hamer - pdh@nilu.no
!*
!* Unless explicitly acquired and licensed from Licensor under another license,
!* the contents of this file are subject to the Reciprocal Public License ("RPL")
!* Version 1.5, https://opensource.org/licenses/RPL-1.5 or subsequent versions as
!* allowed by the RPL, and You may not copy or use this file in either source code
!* or executable form, except in compliance with the terms and conditions of the RPL.
!*
!* All software distributed under the RPL is provided strictly on an "AS IS" basis,
!* WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, AND LICENSOR HEREBY
!* DISCLAIMS ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF
!* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT, OR NON-INFRINGEMENT.
!* See the RPL for specific language governing rights and limitations under the RPL.
!*
!* ==========================================
!* The dispersion model EPISODE (Grønskei et. al., 1993; Larssen et al., 1994;
!* Walker et al., 1992, 1999; Slørdal et al., 2003, 2008) is an Eulerian grid model
!* with embedded subgrid models for calculations of pollutant concentrations resulting
!* from different types of sources (area-, line- and point sources). EPISODE solves the
!* time dependent advection/-diffusion equation on a 3 dimensional grid.
!* Finite difference numerical methods are applied to integrate the solution forward in time.
!* It also includes extensions as the implementation of a simplified EMEP photochemistry
!* scheme for urban areas (Walker et al. 2004) and a scheme for Secondary Organic Aerosol
!* implemented by Håvard Slørdal
!*
!* Grønskei, K.E., Walker, S.E., Gram, F. (1993) Evaluation of a model for hourly spatial
!*    concentrations distributions. Atmos. Environ., 27B, 105-120.
!* Larssen, S., Grønskei, K.E., Gram, F., Hagen, L.O., Walker, S.E. (1994) Verification of
!*    urban scale time-dependent dispersion model with sub-grid elements in Oslo, Norway.
!*    In: Air poll. modelling and its appl. X. New York, Plenum Press.
!* Slørdal, L.H., McInnes, H., Krognes, T. (2008): The Air Quality Information System AirQUIS.
!*    Info. Techn. Environ. Eng., 1, 40-47, 2008.
!* Slørdal, L.H., Walker, S.-E., Solberg, S. (2003) The Urban Air Dispersion Model EPISODE
!*    applied in AirQUIS. Technical Description. NILU TR 12/2003. ISBN 82-425-1522-0.
!* Walker, S.E., Grønskei, K.E. (1992) Spredningsberegninger for on-line overvåking i Grenland.
!*    Programbeskrivelse og brukerveiledning. Lillestrøm,
!*    Norwegian Institute for Air Research (NILU OR 55/92).
!* Walker, S.E., Slørdal, L.H., Guerreiro, C., Gram, F., Grønskei, K.E. (1999) Air pollution
!*    exposure monitoring and estimation. Part II. Model evaluation and population exposure.
!*    J. Environ. Monit, 1, 321-326.
!* Walker, S.-E., Solberg, S., Denby, B. (2003) Development and implementation of a
!*    simplified EMEP photochemistry scheme for urban areas in EPISODE. NILU TR 13/2013.
!*    ISBN 82-425-1524-7
!*****************************************************************************!

      subroutine CPHOTR(ICV)

! The subroutine calculates photochemistry for all irregular
! receptor points.
! ----------------------------------------------------------------------------------
! Based on:
! Version Episode 5.5 (May29, 2012) prepared for BB-Stand-Alone
! Original source code of EPISODE by Sam-Erik Walker (NILU)
!
! Sam-Erik Walker
! Norwegian Institute for Air Research (NILU)
! Instituttveien 18 P.O. Box 100
! N-2027 Kjeller, NORWAY
! Tel: +47 63898000 Fax: +47 63898050
! E-mail: sam-erik.walker@nilu.no
!
! ----------------------------------------------------------------------------------
! ----------------------------------------------------------------------------------
! REVISION HISTORY
!
!           2016  M. Karl: modified setting of cloud cover.
!                          Lower concentration cut off for very small values.
!                          Modified call to photo.
!                          Added debugging statements for testing
!
! ----------------------------------------------------------------------------------

      use mod_util
      use mod_main
      use mod_site
      use mod_time
      use mod_mete
      use mod_conc
      use mod_depo
      use mod_phot

      implicit none

      integer :: ICV

! ICV - Update concentration indicator

! Local variables

      real :: CLOUV
      real :: CV(3)
      real :: LATIV
      real :: TAIRV
      real :: XRV
      real :: YRV
      real :: ZRV
      integer :: DAYMV
      integer :: HOURV
      integer :: IR
      integer :: IX
      integer :: IY
      integer :: IZ
      integer :: MNTHV
      integer :: YEARV

!MSK start
      integer :: ic
      integer :: SECOV
      integer :: MINUV
      real    :: SunHr
      integer :: LTOGV
      integer :: DAYYV

!         real   :: DecAng
!         real   :: CosZen
!         real   :: ZenAng
!         real   :: SinSH
!         real, parameter :: pii   = 3.1415927
!         real, parameter :: radd  = pii/180.0
!MSK end

! CLOUV - Cloud cover value
! CV    - Receptor concentration values
! LATIV - Latitude value
! TAIRV - Air temperature value
! XRV   - Receptor point x-coordinate
! YRV   - Receptor point y-coordinate
! ZRV   - Receptor point z-coordinate
! DAYMV - Day value
! HOURV - Hour value
! IR    - Receptor point index
! IX    - Main grid index in x-direction
! IY    - Main grid index in y-direction
! IZ    - Main grid index in z-direction
! MNTHV - Month value
! YEARV - Year value

! Check photochemistry indicator

!MSK      IF (PHOTPS .NE. 1) RETURN

! Initialize latitude and date/time

      LATIV = SITELA

      YEARV = YEAR
      MNTHV = MNTH
      DAYMV = DAYM
      HOURV = HOUR

!MSK start
      DAYYV = DAYY
      MINUV = MINU
      SECOV = SECO
      LTOGV = LTOGMT
!MSK end

! Go through all receptor points

      DO 100 IR = 1,NR

! Define receptor point coordinates

          XRV = XR(IR)
          YRV = YR(IR)
          ZRV = ZR(IR)

! Get main grid indices

          call getmgi(1,XRV,YRV,ZRV,IX,IY,IZ)

! Set cloud cover

!MSK start
! Estimation of cloud cover:
! If total solar radiation is available, then calculate
! else take from cloud cover file (constant value)

          IF (TSRADFE) THEN
             if (TSRAD(IX,IY).eq.0.0) then
                ! Nighttime cloud cover not relevant
                CLOUV = 1.0
            !    ZenAng = 1.0
            !    SinSH  = 0.0
             else
                SunHr = HOURV + (60.*MINUV + SECOV)/3600.

            ! ! Calculate the Zenith angle for test purpose
            ! DecAng = 23.27*COS(2.0*pii*(DAYMV - 172.0)/365.0)
            ! CosZen = COS(radd*LATIV)*COS(radd*DecAng) *          &
            !         COS(2.0*pii*(SunHr - 12.0)/24.0) +           &
            !         SIN(radd*LATIV)*SIN(radd*DecAng)
            ! ZenAng = ACOS(max(0.017,CosZen))/radd
            ! SinSH  = SIN(radd*ZenAng)

                CLOUV = cloudfraction(SunHr,DAYMV,LATIV,TSRAD(IX,IY))
             endif
          ELSE
             CLOUV = CLOU(IX,IY)
          ENDIF
         ! if (ir.le.30) print *,'cphotr CLOUV:',ir,TSRAD(IX,IY),SinSH,CLOUV
!MSK end

! Set air temperature

          TAIRV = TAIR(IX,IY) + CTOK

!MSK start
! Concentration cut off for very small
          do IC=1,NC
             if ( CR(IC,IR) .lt. 1.e-20 )   CR(IC,IR) = 0.0
          enddo
!MSK end

!MSK start
!MSK: for receptor points and lines only photochemical balance
!
! Mapping photochemistry indices to compound indices

!MSK debug start
!MSK debug:   Check NO and NO2 before photo
!MSK debug          if (ir.eq.10000) then
!MSK debug             print *,'cphotr NO  before photo',CR(2,ir)
!MSK debug             print *,'cphotr NO2 before photo',CR(3,ir)
!MSK debug             print *,'cphotr PM2.5 before photo',CR(15,ir)
!MSK debug          endif
!MSK debug end

! PHOTPS = 1 :
! 1 = NO2
! 2 = NO
! 3 = O3

! Copy NO2, NO and O3 concentrations (ug/m3)

          do IC=1,NC
            if(trim(CMPND(IC)) == 'NO2')   CV( 1 )  = CR(IC,IR)
            if(trim(CMPND(IC)) == 'NO')    CV( 2 )  = CR(IC,IR)
            if(trim(CMPND(IC)) == 'O3')    CV( 3 )  = CR(IC,IR)
          enddo
!MSK end.

! Calculate photochemical balance between NO2, NO and O3

!MSK          CALL PHOTO(LATIV,YEARV,MNTHV,DAYMV,HOURV,CLOUV,TAIRV,CV)
            call photo(LATIV,LTOGV,YEARV,MNTHV,DAYMV,DAYYV,HOURV,MINUV,SECOV,CLOUV,TAIRV,CV)

          !  print *,'cphotr after photo ',CLOUV,TAIRV,CV

! Copy NO2, NO and O3 concentrations back (ug/m3)

!MSK start

          do IC=1,NC
            if(trim(CMPND(IC)) == 'NO2')   CR(IC,IR) = CV( 1 )
            if(trim(CMPND(IC)) == 'NO')    CR(IC,IR) = CV( 2 )
            if(trim(CMPND(IC)) == 'O3')    CR(IC,IR) = CV( 3 )
          enddo

!MSK end

!MSK debug start
!MSK debug:   Check NO and NO2 after photo
!MSK debug          if (ir.eq.10000) then
!MSK debug             print *,'cphotr NO  after photo',CR(2,ir)
!MSK debug             print *,'cphotr NO2 after photo',CR(3,ir)
!MSK debug             print *,'cphotr PM2.5 before photo',CR(15,ir)
!MSK debug          endif
!MSK debug end

! Next receptor point

  100 CONTINUE

      RETURN

! End of subroutine CPHOTR

      end subroutine cphotr

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


text:
C--------------------------------------------------------------------------
        subroutine ftgcno(iunit,casesn,templt,colnum,status)

C       determine the column number corresponding to an input column name.
C       This supports the * and ? wild cards in the input template.

C       iunit   i  Fortran i/o unit number
C       casesn  l  true if an exact case match of the names is required
C       templt  c  name of column as specified in a TTYPE keyword
C       colnum  i  number of the column (first column = 1)
C                       (a value of 0 is returned if the column is not found)
C       status  i  returned error status

C       modified by Wm Pence, HEASARC/GSFC, December 1994

        integer iunit,colnum,status
        character*(*) templt
        logical casesn
        character*8 dummy

        call ftgcnn(iunit,casesn,templt,dummy,colnum,status)
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine usersections(inpc,textpart,set,istartset,iendset,
     &  ialset,nset,ielmat,matname,nmat,
     &  irstrt,istep,istat,n,iline,ipol,
     &  inl,ipoinp,inp,lakon,ielprop,nprop,nprop_,prop,
     &  ipoinpc,mi,ier)
!
!     reading the input deck: *USER SECTION
!
      implicit none
!
      character*1 inpc(*)
      character*8 lakon(*)
      character*80 matname(*),material
      character*81 set(*),elset
      character*132 textpart(16)
!
      integer mi(*),istartset(*),iendset(*),ialset(*),id,
     &     ielmat(mi(3),*),irstrt(*),nset,nmat,ndprop,npropstart,
     &     istep,istat,n,key,i,j,k,imaterial,ipos,lprop,ipoinpc(0:*),
     &     iline,ipol,inl,ipoinp(2,*),inp(3,*),ielprop(*),nprop,nprop_,
     &     iset,ier
!
      real*8 prop(*)
!
      if((istep.gt.0).and.(irstrt(1).ge.0)) then
        write(*,*)
     &       '*ERROR reading *USER SECTION: *USER SECTION should'
        write(*,*) '  be placed before all step definitions'
        ier=1
        return
      endif
!
      do i=2,n
        if(textpart(i)(1:9).eq.'MATERIAL=') then
          material=textpart(i)(10:89)
        elseif(textpart(i)(1:6).eq.'ELSET=') then
          elset=textpart(i)(7:86)
          elset(81:81)=' '
          ipos=index(elset,' ')
          elset(ipos:ipos)='E'
!
        elseif(textpart(i)(1:10).eq.'CONSTANTS=') then
          read(textpart(i)(11:20),'(i10)',iostat=istat) ndprop
          if(istat.gt.0) then
            call inputerror(inpc,ipoinpc,iline,
     &           "*USER SECTION%",ier)
            return
          endif
        else
          write(*,*)
     &       '*WARNING reading *USER SECTION: parameter not recognized:'
          write(*,*) '         ',
     &         textpart(i)(1:index(textpart(i),' ')-1)
          call inputwarning(inpc,ipoinpc,iline,
     &         "*USER SECTION%")
        endif
      enddo
!
!     check for the existence of the set and the material
!
      do i=1,nmat
        if(matname(i).eq.material) exit
      enddo
      if(i.gt.nmat) then
        write(*,*)
     &       '*ERROR reading *USER SECTION: nonexistent material'
        write(*,*) '  '
        call inputerror(inpc,ipoinpc,iline,
     &       "*USER SECTION%",ier)
        return
      endif
      imaterial=i
!
c      do i=1,nset
c        if(set(i).eq.elset) exit
c      enddo
      call cident81(set,elset,nset,id)
      i=nset+1
      if(id.gt.0) then
        if(elset.eq.set(id)) then
          i=id
        endif
      endif
      if(i.gt.nset) then
        elset(ipos:ipos)=' '
        write(*,*) '*ERROR reading *USER SECTION: element set ',elset
        write(*,*) '  has not yet been defined. '
        call inputerror(inpc,ipoinpc,iline,
     &       "*USER SECTION%",ier)
        return
      endif
      iset=i
!
      npropstart=nprop
!
!     reading the element properties
!
      if(ndprop.gt.0) then
!
!     general case
!
        lprop=0
!
        do
          call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &         ipoinp,inp,ipoinpc)
          if((istat.lt.0).or.(key.eq.1)) exit
          do k=1,n
            lprop=lprop+1
            if(lprop.gt.ndprop) exit
            read(textpart(k),'(f20.0)',iostat=istat)
     &           prop(nprop+lprop)
            if(istat.gt.0) then
              call inputerror(inpc,ipoinpc,iline,
     &             "*USER SECTION%",ier)
              return
            endif
          enddo
        enddo
        nprop=nprop+ndprop
      else
        call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &       ipoinp,inp,ipoinpc)
      endif
!
      if(nprop.gt.nprop_) then
        write(*,*) '*ERROR reading *USER SECTION: increase nprop_'
        ier=1
        return
      endif
!
!     assigning the elements of the set the appropriate material
!     and property pointer
!
      do j=istartset(iset),iendset(iset)
        if(ialset(j).gt.0) then
          if(lakon(ialset(j))(1:1).ne.'U') then
            write(*,*) '*ERROR reading *USER SECTION: element ',
     &           ialset(j),' is no user element.'
            ier=1
            return
          endif
!
          ielmat(1,ialset(j))=imaterial
          if(ndprop.gt.0) ielprop(ialset(j))=npropstart
        else
          k=ialset(j-2)
          do
            k=k-ialset(j)
            if(k.ge.ialset(j-1)) exit
            if(lakon(k)(1:1).ne.'U') then
              write(*,*) '*ERROR reading *USER SECTION: element ',
     &             k,' is no user element.'
              ier=1
              return
            endif
            ielmat(1,k)=imaterial
            if(ndprop.gt.0) ielprop(k)=npropstart
          enddo
        endif
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text: junior at UC Davis UCD, no job, just chilling

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


program cobldrv
    use param
    use complex_oblate_swf

    real(knd)       x, ca, arg1, darg, step1, step2, step3, xneu
    complex(knd)    c
    character       chr

    integer         i, im, j, mmin, minc, mnum, m, l, lnum, ioprad, iopang, &
                    iopnorm, ioparg, narg, kind, kindd, kindq

    real (knd), dimension(:), allocatable ::        eta
    complex(knd), dimension(:), allocatable ::      r1c, r1dc, r2c, r2dc
    integer, dimension(:), allocatable ::           ir1e, ir1de, ir2e, ir2de, naccr
    complex(knd), dimension (:,:), allocatable ::   s1c, s1dc
    integer, dimension(:,:), allocatable ::         is1e, is1de, naccs, naccds

    kindd =  8
    kindq = 16

    open(1, file='coblfcn.dat')
    open(20, file='fort.20')
    open(30, file='fort.30')
    open(40, file='fort.20')
    open(50, file='fort.30')
    open(60, file='fort.20')

    read(1,*) mmin, minc, mnum, lnum
        read(1,*) ioprad, iopang, iopnorm
        read(1,*) c, x
        if(iopang.ne.0) read(1,*) ioparg, arg1, darg, narg

    allocate (eta(narg), r1c(lnum), r1dc(lnum), r2c(lnum), r2dc(lnum))
    allocate (ir1e(lnum), ir1de(lnum), ir2e(lnum), ir2de(lnum), naccr(lnum))
    allocate (s1c(lnum, narg), s1dc(lnum, narg))
    allocate (is1e(lnum, narg), is1de(lnum, narg), naccs(lnum, narg), naccds(lnum, narg))

    if(knd.eq.kindd) minacc = 8
    if(knd == kindq .and. aimag(c) <= 20.0e0_knd) minacc = 15
    if(knd == kindq .and. aimag(c) >  20.0e0_knd) minacc =  8

    ca=abs(c)
    ndec=precision(ca)
    nex=range(ca) - 1

    nbp = int(2 * (real(c) + aimag(c)) / 3.1416)
    if(lnum < nbp .and. (2 * (lnum / 2) /= lnum)) lnum = lnum + 1
    maxe = max(50, nbp + 30) + 10
    maxe2 = maxe + maxe
    if(maxe2 < lnum) maxe2 = lnum
    maxm = mmin + minc * (mnum - 1)
    maxlp = lnum + maxm + 1
    maxint = 2 * (nbp + 33) + 306
    maxj = lnum + 3 * ndec + int(ca) + 5 + maxm
    maxp = max(lnum + 3 * ndec + int(ca) + 5, maxlp + 5)
    maxn = maxj
    maxpdr = 2 * int(ca) + 4 * ndec + int(100 * x) + 8
    neta = 30
    step1 = 0.1e0_knd
    step2 = 0.1e0_knd
    step3 = 0.8e0_knd
    nstep1 = 1
    nstep2 = 1
    nstep3 = 3
    ngau = 100 * (2 + int(ca / 500.0e0_knd))
    xneu = 0.3e0_knd

    if(ioprad == 2) then
        if(x < 0.2e0_knd) then
            step1 = x / 4.0e0_knd
            step2 = 0.075e0_knd
            step3 = 1.0e0_knd - step1 - step2
            nstep1 = 1
            nstep2 = 1
            nstep3 = 4
            ngau=100 * (2 + int(ca / 500.0e0_knd))
        end if

        if(x.lt.0.1e0_knd) then
            step1 = x /4.0e0_knd
            step2 = 0.025e0_knd
            step3 = 1.0e0_knd - step1 - step2
            nstep1 = 1
            nstep2 = 1
            nstep3 = 4
            ngau = 200
            if(ca > 500.0e0_knd) ngau = 200 * (2 + int(ca /1000.0e0_knd))
        end if

        if(x < 0.05e0_knd) then
            step1 = x / 15.0e0_knd
            step2 = 0.002e0_knd
            step3 = 1.0e0_knd - step1 - step2
            nstep1 = 1
            nstep2 = 1
            nstep3 = 2
            ngau = 300
            if(ca >  500.0e0_knd) ngau =  500
            if(ca > 1000.0e0_knd) ngau =  800
            if(ca > 1500.0e0_knd) ngau = 1000
            if(ca > 2000.0e0_knd) ngau = 1200
            if(ca > 2500.0e0_knd) ngau = 1500
            if(ca > 3000.0e0_knd) ngau = 1700
            if(ca > 3500.0e0_knd) ngau = 1900
            if(ca > 4000.0e0_knd) ngau = 2200
            if(ca > 4500.0e0_knd) ngau = 2500
            if(ca > 5000.0e0_knd) ngau = 2500 + 300 * int((ca - 4500.0e0_knd) / 500.0e0_knd)
        end if

        if(x < 0.01e0_knd) then
            step1 = x / 15.0e0_knd
            step2 = 0.002e0_knd
            step3 = 1.0e0_knd - step1 - step2
            nstep1 = 1
            nstep2 = 1
            nstep3 = 2
            ngau = 300
            if(ca >  500.0e0_knd) ngau =  600
            if(ca > 1000.0e0_knd) ngau =  800
            if(ca > 1500.0e0_knd) ngau = 1000
            if(ca > 2000.0e0_knd) ngau =  300 * int(ca /500.0e0_knd)
        end if

        if(x < 0.001e0_knd) then
            step1 = x / 15.0e0_knd
            step2 = 0.002e0_knd
            step3 = 1.0e0_knd - step1 - step2
            nstep1 = 3
            nstep2 = 1
            nstep3 = 2
            ngau = 600
            if(ca >  300.0e0_knd) ngau =  800
            if(ca > 1000.0e0_knd) ngau =  900
            if(ca > 1500.0e0_knd) ngau = 1000
            if(ca > 2000.0e0_knd) ngau = 300 * int(ca / 500.0e0_knd)
            if(aimag(c) > 5.0e0_knd) ngau = ngau + (ngau / 5) * min(5, int(aimag(c)) - 5)
            if(knd.eq.kindq.and.x.lt.0.00001e0_knd) ngau=2*ngau
            if(knd.eq.kindq.and.x.lt.0.000001e0_knd) ngau=2*ngau
            if(knd.eq.kindd.and.aimag(c).gt.4.0e0_knd.and.x.lt. 0.00001e0_knd) ngau=2*ngau
            if(knd.eq.kindd.and.aimag(c).gt.4.0e0_knd.and.x.lt. 0.000001e0_knd) ngau=2*ngau
        end if

        xneu = 0.3e0_knd
        if(ca > 100.0e0_knd) xneu = 0.04e0_knd
        if(ca > 600.0e0_knd) xneu = 0.03e0_knd
        if(ca > 800.0e0_knd) xneu = 0.01e0_knd
        if(aimag(c) > 50.0e0_knd) xneu = 0.01e0_knd

        if(x >= 0.01e0_knd) then
            if(knd == kindd) then
                if(x >= 0.01e0_knd) maxn = 2 * int(25 / (x * x) + 300 / x + 3 * ca + 1250 * knd) + 5
                if(x >= 0.1e0_knd)  maxn = 2 * int((lnum + ca / 5 + 0.5e0_knd * maxm + 200) * 1.4e0_knd / x) + 5
                if(x >= 0.5e0_knd)  maxn = 2 * int((lnum + ca / 5 + 0.5e0_knd * maxm + 300) / x) + 5
                if(x >= 1.0e0_knd)  maxn = 2 * int(lnum + ca / 5 + 0.5e0_knd * maxm + 300) + 5
            end if
            if(knd.eq.kindq) then
                if(x >= 0.01e0_knd) maxn = 2 * int(25 / (x * x) + 400 / x + 3 * ca + 1250 * knd) + 5
                if(x >= 0.1e0_knd)  maxn = 2 * int((lnum + ca / 5 + 0.5e0_knd * maxm+ 350) * 1.4e0_knd / x) + 5
                if(x >= 0.5e0_knd)  maxn = 2 * int((lnum + ca / 5 + 0.5e0_knd * maxm + 400) / x) + 5
                if(x >= 1.0e0_knd)  maxn = 2 * int(lnum + ca / 5 + 0.5e0_knd * maxm + 400) + 5
                end if
            maxn = maxn + maxm
        end if
    end if

    maxp = max(maxn, maxp, maxpdr)
    maxq = lnum + 3 * ndec + int(ca) + maxm + maxm + 4
    if(knd == kindd .and. aimag(c) < 10.0e0_knd .and. real(c) <= 60.0e0_knd &
        .and. real(c) >= 10.0e0_knd .and. mmin <= 40 .and. x <= 0.99e0_knd &
        .and. x > 0.1e0_knd) maxq = max(maxq, 250 - int(50 * x) + 4 + maxm + maxm)
    maxdr = maxpdr / 2 + 1
    maxd = maxn / 2 + 1
    maxmp = maxm + maxm + 5
    maxt = 1
    jnebmax = 30
    jnenmax = 10
    if(x < 0.05e0_knd) jnenmax = 1
    if(iopang /= 0) maxt = narg

    if (iopang /= 0) then
        do j = 1, narg
            eta(j) = arg1 + (j - 1) * darg
        end do
    end if

    do im = 1, mnum
        m = mmin + (im - 1) * minc

        call coblfcn(c, m, lnum, ioprad, x, iopang, iopnorm, narg, eta, &
                     r1c, ir1e, r1dc, ir1de, r2c, ir2e, r2dc, ir2de, naccr, &
                     s1c, is1e, s1dc, is1de, naccs, naccds)

        if (ioprad /= 0) then

            if(knd == kindd) write(20, 10) x, c, m
10          format(1x,e23.14,e23.14,e23.14,i5)
            if(knd == kindq) write(20,20) x, c, m
20          format(1x,e39.30,e39.30,e39.30,i5)

            do i = 1, lnum
                l = m + i - 1
                if (ioprad == 1) write(20, 30) l, r1c(i), ir1e(i), r1dc(i), ir1de(i), naccr1, ' '

                if (ioprad == 2) then
                    chr = 'w'
                    if (naccr(i) < 0) chr = 'e'
                    write(20, 40) l, r1c(i), ir1e(i), r1dc(i), ir1de(i), r2c(i), ir2e(i), r2dc(i), ir2de(i), abs(naccr(i)), chr
                end if

30              format(1x, i5, 2x, 2(f17.14, 1x, f17.14, i6, 2x), i2, a)
40              format(1x, i5, 2x, 2(f17.14, 1x, f17.14, i6, 2x), /, 8x, 2(f17.14, 1x, f17.14, i6, 2x), i2, a)
            end do
        end if

        if (iopang /= 0) then

            if(knd == kindd) write(30, 50) c, m
            if(knd == kindq) write(30, 60) c, m
50          format(1x,'c = ',e23.14, e23.14,'; m = ',i5)
60          format(1x,'c = ',e39.30, e39.30,'; m = ',i5)

            do i = 1, lnum
                l = m + i - 1

                write(30, "(1x, i6)") l

                do j = 1, narg
                    arg = arg1 + (j - 1) * darg

                    if(ioparg.eq.1.and.iopang.eq.1) write(30, 70) eta(j), s1c(i,j), is1e(i,j), naccs(i,j)
                    if(ioparg.eq.1.and.iopang.eq.2) write(30, 80) eta(j), s1c(i,j), is1e(i,j), s1dc(i,j), is1de(i,j), naccs(i,j), naccds(i,j)

70                  format(1x, f19.14, 2x, f17.14, 1x, f17.14, 2x, i5, 2x, ', ', i2)
80                  format(1x, f19.14, 2x, f17.14, 1x, f17.14, 2x, i5, 2x, f17.14, 1x, f17.14, 2x, i5, 2x, i2, ', ', i2)

                end do
            end do
        end if
    end do

    deallocate (is1e, is1de, naccs, naccds)
    deallocate (s1c, s1dc)
    deallocate (ir1e, ir1de, ir2e, ir2de, naccr)
    deallocate (eta, r1c, r1dc, r2c, r2dc)
    close(60)
    close(50)
    close(40)
    close(30)
    close(20)
    close(1)

end program cobldrv


text:
http://www.upte.org/ University Professional and Technical Employees (CWA 9119) is the Labor Unions union of technical and professional employees at the University of California. We are Staff Research Associates, Computer Resource Specialists, Clinical Lab Techs, Editors, Student Affairs Officers, Social Workers, Writers, Museum Scientists and Lab Assistants, among may other titles. There are more than 11,000 UC employees covered by UPTE contracts. UPTE was founded in 1990 by a group of employees who believed that UC workers would benefit by having a union to safeguard and expand our rights in the workplace. There are UPTE locals at all UC campuses, medical centers, and the national labs administered by UC. In 1993, UPTE members voted to affiliate with the http://www.cwaunion.org/ Communications Workers of America (CWA), a 700,000 member union in the http://www.aflcio.org/ AFLCIO, in order to better organize and represent UC workers.

UPTE works to protect and improve the rights of employees. We are willing to stand up to the University of California to ensure that professional and technical employees receive the respect and fairness they deserve. UPTE represents employees who face problems on the job, including unfair layoffs and performance evaluations, disciplinary actions, and much more. UPTE members in the technical, research professional and health care professional units bargain with UC over wages and working conditions.

See also: UPTE Strike 2009


!---------------------------------------------------!
! Copyright (c) 2017 Shunsuke A. Sato               !
! Released under the MIT license                    !
! https://opensource.org/licenses/mit-license.php   !
!---------------------------------------------------!
!-------10--------20--------30--------40--------50--------60--------70--------80--------90
program main
  use global_variables
  implicit none
  integer :: it,ikz,ikr
  real(8) :: jz_intra,jz_inter,Etz

  call MPI_init(ierr)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,Nprocs,ierr)
  call MPI_COMM_RANK(MPI_COMM_WORLD,Myrank,ierr)

  call input
  call preparation
  call init_Ac

  call time_propagation

  if(Myrank == 0) write(*,*)  'This calculation is shutdown successfully!'
  call MPI_FINALIZE(ierr)
end program main


text:
C    @(#)out_bus.f	20.4 2/13/96
C****************************************************************
C
C   File: out_bus.f
C   Purpose: Routine to obtain BCD image of INPUT/OUTPUT/SYSTEM
C            data
C
C   Author: Walt Powell  Date: 20 February 1992
C                        Modified: 20 February 1992
C   Called by:
C
C****************************************************************
C
        subroutine out_bus (nb, datarec)

        character datarec * (*)
c
c       This subroutine returns WSCC-formated output data records.
c       Output parameter:
c
c       datarec - a character string for storing data
c
        include 'ipfinc/parametr.inc'

        include 'ipfinc/alpha.inc'
        include 'ipfinc/arcntl.inc'
        include 'ipfinc/area.inc'
        include 'ipfinc/blank.inc'
        include 'ipfinc/branch.inc'
        include 'ipfinc/bus.inc'
        include 'ipfinc/cbus.inc'
        include 'ipfinc/lfiles.inc'
        include 'ipfinc/prt.inc'
        include 'ipfinc/xdata.inc'
        include 'ipfinc/dc2t.inc'
        include 'ipfinc/dcmt.inc'
        include 'ipfinc/ordsta.inc'

        character type * 1

        call bcdbus (nb, datarec)
        if (ordvlt .eq. 1) then
           kt = nb
        else
           kt = inp2opt(nb)
        endif
	ntyp = ntypu(kt)
        call typno (type, ntyp)

	if (ntyp .eq. 5) then
           do idc = 1, kdtot
               if (dc2t(1,idc) .eq. nb ) go to 10
               if (dc2t(3,idc) .eq. nb ) go to 20
           enddo

   10      pdc = dc2t(42,idc)
           qdc = dc2t(44,idc)
           pvlv = abs(pdc - 0.001 * dc2t(40,idc) * dc2t(39,idc))
           qvlv = qdc
           dci = dc2t(39,idc)
           angle = 57.2957795 * dc2t(22,idc)
           vdc = dc2t(40,idc)
           go to 30

   20      pdc = dc2t(43,idc)
           qdc = dc2t(45,idc)
           pvlv = abs(pdc + 0.001 * dc2t(41,idc) * dc2t(39,idc))
           qvlv = qdc
           dci = -dc2t(39,idc)
           angle = 57.2957795 * dc2t(26,idc)
           vdc = dc2t(41,idc)

   30      write( datarec(21:), 40 )  pdc, qdc, vdc, angle, pvlv,
     1                                   qvlv
   40      format(6e15.7)

        else if (ntyp .eq. 12) then

           do idc = 1,mtdcbs
               if (dcmtbs(1,idc) .eq. nb ) go to 110
           enddo
           call erexit
  110      if (dcmtbs(3,idc) .eq. 0) then
              pdc = 0.0
              qdc = 0.0
              pvlv = 0.0
              qvlv = 0.0
              dci = 0.0
	      angle = 0.0
              vdc = 0.0
           else
              pdc = dcmtbs(25,idc)
              qdc = dcmtbs(26,idc)
              pvlv = abs(pdc - dcmtbs(19,idc))
              qvlv = qdc
              if (dcmtbs(20,idc) .ne. 0.0) then
                 dci = 1000.0 * dcmtbs(19,idc) / dcmtbs(20,idc)
              else
                 dci = 0.0
              endif
              angle = 57.2957795 * dcmtbs(13,idc)
              vdc = dcmtbs(24,idc)
           endif
           write( datarec(21:), 40 )  pdc, qdc, vdc, angle, pvlv,
     1                                   qvlv
        else
           call nrpqv (kt, pk, dpk, qk, dqk, vk)

           pload = (ploadu(kt) + inetr(kt) * vk) * bmva
           qload = (qloadu(kt) - ineti(kt) * vk) * bmva
           pgen = (pk - inetr(kt) * vk) * bmva + pload
           voltkv = vk * base(nb)
           if (e(kt) .eq. 0.0 .and. f(kt) .eq. 0.0) then
              degree = 0.0
           else
              degree = 57.2957795 * datan2( f(kt), e(kt) )
           endif
           pvlv = 0.0
           qvlv = 0.0
           call typno( type, ntyp )
           call allocq (nb, qk, qgen, qgnmax, qgnmin, qld, totcap,
     1                  usecap, totrek, userek, unsked, qerr)
           if ( usecap + unsked .gt. 0.0 ) then
              if ( totcap .gt. 0.0 ) then
                 bsched = totcap
              else
                 bsched = totrek
              end if
              bused = usecap + unsked
           else if ( userek + unsked .lt. 0.0 ) then
              if ( totrek .lt. 0.0 ) then
                 bsched = totrek
              else
                 bsched = totcap
              end if
              bused = userek + unsked
           else if ( unsked .lt. 0.0 ) then
              if ( totrek .lt. 0.0 ) then
                 bsched = totrek
              else
                 bsched = totcap
              end if
              bused = unsked
           else if ( unsked .gt. 0.0 ) then
              if ( totcap .gt. 0.0 ) then
                 bsched = totcap
              else
                 bsched = totrek
              end if
              bused = unsked
           else if ( totrek .lt. 0.0 .and. totcap .eq. 0.0 ) then
              bsched = totrek
              bused = userek
           else
              bsched = totcap
              bused = usecap
           end if
           write( datarec(21:), 160 )  pgen, qgen, voltkv, degree,
     &                                 pload, qload, bused, bsched,
     &                                 usecap, totcap, userek, totrek,
     &                                 unsked
  160      format(13e15.7)
	endif
	return
	end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      SUBROUTINE CSTAB2(ITERM)
*
*
*       Degenerate triple chain stability test.
*       ---------------------------------------
*
      INCLUDE 'commonc.h'
      INCLUDE 'common2.h'
      COMMON/CHREG/  TIMEC,TMAX,RMAXC,CM(10),NAMEC(6),NSTEP1,KZ27,KZ30
      REAL*8  M,MB,MB1,R2(NMX,NMX),XCM(3),VCM(3),XX(3,3),VV(3,3),
     &        XCM3(3),VCM3(3)
      INTEGER  IJ(NMX)
*
*
*       Sort particle separations (I1 & I2 form closest pair).
      CALL R2SORT(IJ,R2)
      I1 = IJ(1)
      I2 = IJ(2)
      I3 = IJ(3)
      I4 = IJ(4)
      MB = M(I1) + M(I2)
      MB1 = MB + M(I3)
      MB2 = MB1 + M(I4)
*
*       Form output diagnostics with smallest binary as one body.
      VREL2 = 0.0D0
      VREL21 = 0.0D0
      RDOT = 0.0D0
      RDOT3 = 0.0D0
      RB0 = 0.0
      RI2 = 0.0D0
      R4 = 0.0D0
      V4 = 0.0
      RDOT4 = 0.0D0
      DO 10 K = 1,3
          J1 = 3*(I1-1) + K
          J2 = 3*(I2-1) + K
          J3 = 3*(I3-1) + K
          J4 = 3*(I4-1) + K
          RB0 = RB0 + (X(J1) - X(J2))**2
          VREL2 = VREL2 + (V(J1) - V(J2))**2
          RDOT = RDOT + (X(J1) - X(J2))*(V(J1) - V(J2))
          XCM(K) = (M(I1)*X(J1) + M(I2)*X(J2))/MB
          VCM(K) = (M(I1)*V(J1) + M(I2)*V(J2))/MB
          RI2 = RI2 + (X(J3) - XCM(K))**2
          VREL21 = VREL21 + (V(J3) - VCM(K))**2
          RDOT3 = RDOT3 + (X(J3) - XCM(K))*(V(J3) - VCM(K))
          XCM3(K) = (M(I3)*X(J3) + MB*XCM(K))/MB1
          VCM3(K) = (M(I3)*V(J3) + MB*VCM(K))/MB1
          R4 = R4 + (X(J4) - XCM3(K))**2
          V4 = V4 + (V(J4) - VCM3(K))**2
          RDOT4 = RDOT4 + (X(J4) - XCM3(K))*(V(J4) - VCM3(K))
          XX(K,1) = XCM(K)
          XX(K,2) = XCM3(K)
          XX(K,3) = X(J4) - XCM3(K)
          VV(K,1) = VCM(K)
          VV(K,2) = VCM3(K)
          VV(K,3) = V(J4) - VCM3(K)
   10 CONTINUE
*
*       Evaluate orbital elements for inner and outer motion.
      RB = SQRT(RI2)
      R4 = SQRT(R4)
      SEMI = 2.0D0/RB - VREL21/MB1
      SEMI = 1.0/SEMI
      ECC = SQRT((1.0D0 - RB/SEMI)**2 + RDOT3**2/(SEMI*MB1))
      SEMI1 = 2.0/R4 - V4/MB2
      SEMI1 = 1.0/SEMI1
      ECC1 = SQRT((1.0D0 - R4/SEMI1)**2 + RDOT4**2/(SEMI1*MB2))
*
*       Skip if innermost triple is not stable (including ECC > 1).
      RB0 = SQRT(RB0)
      SEMI0 = 2.0/RB0 - VREL2/MB
      SEMI0 = 1.0/SEMI0
      IF (SEMI0.LT.0.0.OR.ECC.GT.1.0) GO TO 40
      ECC0 = SQRT((1.0 - RB0/SEMI0)**2 + RDOT**2/(SEMI0*MB))
      PCRIT0 = stability(M(I1),M(I2),M(I3),ECC0,ECC,0.0D0)*SEMI0
*       Allow 10% extra with MA 1999 (proper test done for outer triple).
      PMIN0 = 1.1*SEMI*(1.0 - ECC)
      IF (PMIN0.LT.PCRIT0) GO TO 40
*
*       Form hierarchical stability ratio (Eggleton & Kiseleva 1995).
*     QL = MB1/M(I4)
*     Q1 = MAX(MB/MB1,MB1/MB)
*     Q3 = QL**0.33333
*     Q13 = Q1**0.33333
*     AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
*
*     EK = AR*SEMI*(1.0D0 + ECC)
      PMIN = SEMI1*(1.0D0 - ECC1)
*
*       Obtain the inclination (in radians).
      CALL INCLIN(XX,VV,XCM3,VCM3,ALPHA)
*
*       Replace the EK criterion by the MA 1999 stability formula.
      PC99 = stability(MB,M(I3),M(I4),ECC,ECC1,ALPHA)*SEMI
*
*       Evaluate the general stability function.
      IF (ECC1.LT.1.0) THEN
          NST = NSTAB(SEMI,SEMI1,ECC,ECC1,ALPHA,MB,M(I3),M(I4))
          IF (NST.EQ.0) THEN
              PCRIT = 0.99*PMIN
          ELSE
              PCRIT = 1.01*PMIN
          END IF
      ELSE
          PCRIT = 1.01*PMIN
      END IF
*
*       Check hierarchical stability condition (SEMI1 > 0 => ECC1 < 1).
      ITERM = 0
      IF (PMIN.GT.PCRIT.AND.SEMI.GT.0.0.AND.SEMI1.GT.0.0) THEN
          ITERM = -1
          WRITE (6,20)  ECC, ECC1, SEMI, SEMI1, PMIN, PCRIT, PC99
   20     FORMAT (' CSTAB2    E =',F6.3,'  E1 =',F6.3,'  A =',1P,E8.1,
     &                     '  A1 =',E8.1,'  PM =',E9.2,'  PC =',E9.2,
     &                     '  PC99 =',E9.2)
          WRITE (6,25)  NAMEC(I1), NAMEC(I2), ECC0, SEMI0, PMIN0, PCRIT0
   25     FORMAT (' INNER TRIPLE    NAM E0 A0 PM0 PC0 ',
     &                              2I6,F7.3,1P,3E10.2)
          RI = SQRT(CM(1)**2 + CM(2)**2 + CM(3)**2)
          EMAX = 0.0
          WRITE (81,30)  TIMEC, RI, NAMEC(I3), ECC, ECC1,
     &                   SEMI, SEMI1, PCRIT/PMIN, 180.*ALPHA/3.14, EMAX
   30     FORMAT (F9.5,F5.1,I6,2F6.3,1P,2E10.2,0P,F5.2,F6.1,F6.3)
          CALL FLUSH(81)
      END IF
*
   40 RETURN
*
      END


! { dg-do compile }

subroutine foo (x)
  integer, pointer, intent (in) :: x
  integer :: i
!$omp parallel private (x)		! { dg-error "INTENT.IN. POINTER" }
!$omp end parallel
!$omp parallel do lastprivate (x)	! { dg-error "INTENT.IN. POINTER" }
  do i = 1, 10
  end do
!$omp simd linear (x)			! { dg-error "INTENT.IN. POINTER" }
  do i = 1, 10
  end do
!$omp single
!$omp end single copyprivate (x)        ! { dg-error "INTENT.IN. POINTER" }
end


SUBROUTINE eta_linear_bottom(H,k,x,int_kx,omega,time,eta,etax,etaxx,etat)
!
USE Precision
USE Constants
IMPLICIT NONE
REAL(KIND=long), INTENT(IN)  :: H,k,x,omega,time,int_kx
REAL(KIND=long), INTENT(OUT) :: eta,etax
REAL(KIND=long), INTENT(OUT) :: etaxx,etat
! REAL(KIND=long), OPTIONAL :: y,etay,etayy !Boundary fitted ?
! Local variables
REAL(KIND=long) :: kx, half_H, sinkx
!
kx=int_kx-omega*time
half_H = H*half
sinkx = sin(kx)
!
eta   =  half_H*cos(kx)
etax  = -half_H*k*sinkx
etat  =  omega*half_H*sinkx
etaxx = -eta*(k)**2

END SUBROUTINE eta_linear_bottom


text:
! ******************************************************************
! ******************************************************************
!
! Este programa tenta colocar todos os itens pequenos dentro da caixa
! maior primeiro, com o objetivo de minimizar o numero de caixas total
!
! ******************************************************************
! ******************************************************************

program algencanma

  use packdat

  implicit none

  ! LOCAL SCALARS
  logical      :: hasAppended, useReduction, useHeuristic
  integer      :: allocerr, i, j, it, container, numberOfSolutions, &
       lastContainer, nAllocTrials, nCntnr
  ! Timing variables
  integer(kind=8) :: tini, tend, ticks_ps
  real(kind=8) :: elapsed, f, fbest, remArea, limItArea, area, tiArea

  ! LOCAL ARRAYS
  character(80)         :: filename, solfname
  character(len=15)     :: strtmp
  real(kind=8), pointer :: x(:),lbest(:)
  integer, allocatable  :: lCntnr(:)

  ! EXTERNAL SUBROUTINES
  external :: myevalf,myevalg,myevalh,myevalc,myevaljac,myevalhc, &
              myevalfc,myevalgjac,myevalgjacp,myevalhl,myevalhlp, &
              alevalfc,alevalgjac

  ! Set problem's data

  filename = 'data.txt'

  solfname = 'solution.csv'

  nAllocTrials = 5

  open(99, file=solfname)

  write(99, *) ''

  close(99)

  useReduction = .true.

  ! Check arguments for Heuristic Mode
  useHeuristic = .true.

  if ( COMMAND_ARGUMENT_COUNT() .EQ. 1 ) then

     CALL GET_COMMAND_ARGUMENT(1, strtmp)

     read(strtmp, *) useHeuristic

  end if

  ! Start timing

  call SYSTEM_CLOCK(tini, ticks_ps)

  ! Load data

  call loadData(filename)

  ! Apply heuristics

  numberOfSolutions = 0

  if ( useReduction ) call reduction(numberOfSolutions, solfname)

  it = numberOfSolutions

  area = 0.0D0

  do i = 1, nItems

     area = area + iWidth(i) * iLength(i)

  end do
  write(*,*) 'Items total area:', area

  allocate(lbest(2 * nTItems + 1), x(2 * nTItems + 1), &
           lCntnr(nContainers), stat=allocerr)

  if ( allocerr .ne. 0 ) then

     write(*,*) 'Allocation errors in main program'

     stop

  end if

  ! Start of main loop

  call getAvContainers(lCntnr, nCntnr)

  lastContainer = 1

  do while ( nItems .gt. 0 )

     tiArea = 0.0D0

     do i = 1, nItems

        tiArea = tiArea + iWidth(i) * iLength(i)

     end do

     do container = lastContainer, nCntnr - 1

        call setCurrContainer(lCntnr(container))

        if ( cWidth * cLength .gt. tiArea ) then

           exit

        end if

     end do

     call setCurrContainer(lCntnr(container))

     write(*, *) 'Selected container:', lCntnr(container)

     filename = 'sol1.asy'

     fbest = 1.0D+20

     write(*,*)
     write(*,*) 'Creating initial allocation'
     write(*,*)

     call sortItems(x)

     call initialpoint(x, i * 123455, .false.)

     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!! !
     ! Deal with only one container !
     ! !!!!!!!!!!!!!!!!!!!!!!!!!!!! !

     ! This step reduces the number of items

     call extractBoxes(x)

     write(filename, FMT="('sol',I0.3,'.asy')") it

     limItArea = 1.0D+20

     hasAppended = .true.

     write(*,*)
     write(*,*) 'Solving allocation problem'
     write(*,*)

     if ( useHeuristic ) call drawSol(x, filename, .false.)

     ! If in 'Heuristic Mode' do not go to solver
     do while ( .not. useHeuristic )

        ! Optimize

        fbest = 1.0D+20

        do i = 1, nAllocTrials

           call initialpoint(x, i * 123456, .false.)

           call solver(x, f, nItems, cWidth, cLength, iWidth, iLength)

           if (f .lt. fbest) then

              fbest = f

              do j = 1, 2 * nItems

                 lbest(j) = x(j)

              end do

              call drawSol(lbest, filename, .false.)

           end if

           if ( fbest .le. 1.0D-08 ) exit

        end do

        if ( .not. hasAppended ) exit

        if ( fbest .le. 1.0D-08 ) then

           remArea = cWidth * cLength

           do i = 1, nItems

              remArea = remArea - iWidth(i) * iLength(i)

           end do

           call appendBox(lbest, remArea, limItArea, hasAppended)

        else

           limItArea = iWidth(nItems) * iLength(nItems)

           call removeAppendedBox()

           remArea = cWidth * cLength

           do i = 1, nItems

              remArea = remArea - iWidth(i) * iLength(i)

           end do

           call appendBox(lbest, remArea, limItArea, hasAppended)

        end if

     end do

     if ( container .ne. nCntnr .and. remainingItems() .gt. 0 ) then

        write(*, *) 'Sobraram ', remainingItems(), 'itens, cancelando container', lCntnr(container)

        lastContainer = container + 1

        call resetPacking()

     else

        ! Found a good solution. Save it, remove the packed items and
        ! rebuild the list of possible containers.

        call saveSol(solfname)

        call removeBoxes()

        call getAvContainers(lCntnr, nCntnr)

        lastContainer = 1

        it = it + 1

     end if

     write(*, *)
     write(*, *) 'Remaining items: ', nItems
     write(*, *)

  end do

  ! Stop timing

  call SYSTEM_CLOCK(tend)

  elapsed = REAL(tend - tini) / REAL(ticks_ps)

  call printStats(elapsed)

  deallocate(x,lbest,stat=allocerr)

  if ( allocerr .ne. 0 ) then

     write(*,*) 'Deallocation error in main program'

     stop

  end if

  stop

end program algencanma

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

subroutine solver(x, f, nItems, cWidth, cLength, iWidth, iLength)

  ! This subroutine sets configuration parameters for the nonlinear
  ! programming solver and calls it.

  implicit none

  ! SCALAR ARGUMENTS
  integer :: nItems
  real(kind=8) :: cWidth, cLength, f

  ! ARRAY ARGUMENTS
  real(kind=8) :: iWidth(nItems), iLength(nItems), x(2 * nItems)

  intent(in   ) :: nItems, cWidth, cLength, iWidth, iLength
  intent(out  ) :: f
  intent(inout) :: x

  ! LOCAL SCALARS
  logical      :: checkder
  integer      :: hnnzmax,inform,jcnnzmax,m,n,nvparam,i
  real(kind=8) :: cnorm,efacc,efstain,eoacc,eostain,epsfeas,epsopt, &
       nlpsupn,snorm

  ! LOCAL ARRAYS
  character(len=15) :: strtmp
  character(len=80) :: specfnm,outputfnm,vparam(10)
  logical           :: coded(11),equatn(0),linear(0)
  real(kind=8)      :: l(2 * nItems),lambda(0),u(2 * nItems)

  ! EXTERNAL SUBROUTINES
  external :: myevalf,myevalg,myevalh,myevalc,myevaljac,myevalhc, &
              myevalfc,myevalgjac,myevalgjacp,myevalhl,myevalhlp, &
              alevalfc,alevalgjac

  ! Number of variables

  n = 2 * nItems

  ! Set lower bounds, upper bounds, and initial guess

  do i = 1, nItems

     l(2 * (i - 1) + 1) = 0.0D0

     u(2 * (i - 1) + 1) = cLength - iLength(i)

     l(2 * i)           = 0.0D0

     u(2 * i)           = cWidth - iWidth(i)

  end do

  ! Constraints

  m = 0

  equatn(1:m) = .false.

  lambda(1:m) = 0.0d0

  linear(1:m) = .true.

  ! Coded subroutines

  coded(1:11) = .false.

  coded(7:8)  = .true.

  ! Upper bounds on the number of sparse-matrices non-null elements

  jcnnzmax = 0

  hnnzmax  = 0

  ! Checking derivatives?

  checkder = .false.

  ! Parameters setting

  epsfeas   = 1.0d-10
  epsopt    = 1.0d-10

  efstain   = sqrt( epsfeas )
  eostain   = epsopt ** 1.5d0

  efacc     = sqrt( epsfeas )
  eoacc     = sqrt( epsopt )

  outputfnm = ''
  specfnm   = 'algencan.dat'

  nvparam = 0

  call algencan(myevalf,myevalg,myevalh,myevalc,myevaljac,myevalhc, &
       alevalfc,alevalgjac,myevalgjacp,myevalhl,myevalhlp,jcnnzmax, &
       hnnzmax,epsfeas,epsopt,efstain,eostain,efacc,eoacc,outputfnm, &
       specfnm,nvparam,vparam,n,x,l,u,m,lambda,equatn,linear,coded, &
       checkder,f,cnorm,snorm,nlpsupn,inform)

end subroutine solver

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

subroutine myevalf(n,x,f,flag)

  implicit none

  ! SCALAR ARGUMENTS
  integer,      intent(in)  :: n
  integer,      intent(out) :: flag
  real(kind=8), intent(out) :: f

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in) :: x(n)

  flag = - 1

end subroutine myevalf

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

subroutine myevalg(n,x,g,flag)

  implicit none

  ! SCALAR ARGUMENTS
  integer,      intent(in)  :: n
  integer,      intent(out) :: flag

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: g(n)

  flag = - 1

end subroutine myevalg

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

subroutine myevalh(n,x,hrow,hcol,hval,hnnz,lim,lmem,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical,      intent(out) :: lmem
  integer,      intent(in)  :: lim,n
  integer,      intent(out) :: flag,hnnz

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: hcol(lim),hrow(lim)
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: hval(lim)

  flag = - 1

  lmem = .false.

end subroutine myevalh

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

subroutine myevalc(n,x,ind,c,flag)

  implicit none

  ! SCALAR ARGUMENTS
  integer,      intent(in)  :: ind,n
  integer,      intent(out) :: flag
  real(kind=8), intent(out) :: c

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in)  :: x(n)

  flag = - 1

end subroutine myevalc

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

subroutine myevaljac(n,x,ind,jcvar,jcval,jcnnz,lim,lmem,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical, intent(out) :: lmem
  integer, intent(in)  :: ind,lim,n
  integer, intent(out) :: flag,jcnnz

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: jcvar(lim)
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: jcval(lim)

  flag = - 1

  lmem = .false.

end subroutine myevaljac

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

subroutine myevalhc(n,x,ind,hcrow,hccol,hcval,hcnnz,lim,lmem,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical, intent(out) :: lmem
  integer, intent(in)  :: ind,lim,n
  integer, intent(out) :: flag,hcnnz

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: hccol(lim),hcrow(lim)
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: hcval(lim)

  flag = - 1

  lmem = .false.

end subroutine myevalhc

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

subroutine myevalfc(n,x,f,m,c,flag)

  use packdat, only: evalover, nItems, iLength

  implicit none

  ! SCALAR ARGUMENTS
  integer,      intent(in)  :: m,n
  integer,      intent(out) :: flag
  real(kind=8), intent(out) :: f

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: c(m)

  ! LOCAL SCALARS
  integer :: i, j

  flag = 0

  f = x(2 * nItems + 1)

  do i = 1, nItems

     do j = i + 1, nItems

        f = f + evalover(x, i, j)

     end do

     c(i) = x(2 * (i - 1) + 1) + iLength(i) - x(2 * nItems + 1)

  end do

end subroutine myevalfc

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

subroutine alevalfc(n,x,f,m,c,flag)

  use packdat, only: evalover, nItems, iLength

  implicit none

  ! SCALAR ARGUMENTS
  integer,      intent(in)  :: m,n
  integer,      intent(out) :: flag
  real(kind=8), intent(out) :: f

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: c(m)

  ! LOCAL SCALARS
  integer :: i, j

  flag = 0

  f = 0.0D0

  do i = 1, nItems

     do j = i + 1, nItems

        f = f + evalover(x, i, j)

     end do

  end do

end subroutine alevalfc

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

subroutine myevalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag)

  use packdat, only : evaldover, nItems

  implicit none

  ! SCALAR ARGUMENTS
  logical,      intent(out) :: lmem
  integer,      intent(in)  :: lim,m,n
  integer,      intent(out) :: flag,jcnnz

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: jcfun(lim),jcvar(lim)
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: g(n),jcval(lim)

  ! LOCAL SCALARS
  integer :: i, j

  flag = 0

  lmem = .false.

  do i = 1, 2 * nItems + 1

     g(i) = 0.0D0

  end do

  jcnnz = 0

  do i = 1, nItems

     do j = i + 1, nItems

        call evaldover(x, i, j, g)

     end do

     jcnnz = jcnnz + 1

     jcfun(jcnnz) = i
     jcvar(jcnnz) = 2 * (i - 1) + 1
     jcval(jcnnz) = 1.0D0

     jcnnz = jcnnz + 1

     jcfun(jcnnz) = i
     jcvar(jcnnz) = 2 * nItems + 1
     jcval(jcnnz) = - 1.0D0

  end do

  g(2 * nItems + 1) = 1.0D0

end subroutine myevalgjac

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

subroutine alevalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag)

  use packdat, only : evaldover, nItems

  implicit none

  ! SCALAR ARGUMENTS
  logical,      intent(out) :: lmem
  integer,      intent(in)  :: lim,m,n
  integer,      intent(out) :: flag,jcnnz

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: jcfun(lim),jcvar(lim)
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: g(n),jcval(lim)

  ! LOCAL SCALARS
  integer :: i, j

  flag = 0

  lmem = .false.

  do i = 1, 2 * nItems

     g(i) = 0.0D0

  end do

  do i = 1, nItems

     do j = i + 1, nItems

        call evaldover(x, i, j, g)

     end do

  end do

  jcnnz = 0

end subroutine alevalgjac

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

subroutine myevalgjacp(n,x,g,m,p,q,work,gotj,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical,   intent(inout) :: gotj
  integer,   intent(in)    :: m,n
  integer,   intent(out)   :: flag
  character, intent(in)    :: work

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in)    :: x(n)
  real(kind=8), intent(inout) :: p(m),q(n)
  real(kind=8), intent(out)   :: g(n)

  flag = - 1

end subroutine myevalgjacp

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

subroutine myevalhl(n,x,m,lambda,sf,sc,hlrow,hlcol,hlval,hlnnz,lim,lmem,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical,      intent(out) :: lmem
  integer,      intent(in)  :: lim,m,n
  integer,      intent(out) :: flag,hlnnz
  real(kind=8), intent(in)  :: sf

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: hlcol(lim),hlrow(lim)
  real(kind=8), intent(in)  :: lambda(m),sc(m),x(n)
  real(kind=8), intent(out) :: hlval(lim)

  flag = - 1

end subroutine myevalhl

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

subroutine myevalhlp(n,x,m,lambda,sf,sc,p,hp,goth,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical,      intent(inout) :: goth
  integer,      intent(in)    :: m,n
  integer,      intent(out)   :: flag
  real(kind=8), intent(in)    :: sf

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in)  :: lambda(m),p(n),sc(m),x(n)
  real(kind=8), intent(out) :: hp(n)

  flag = - 1

end subroutine myevalhlp

max_stars_repo_name: fsobral/pack-u

max_issues_repo_name: fsobral/pack-u

max_forks_repo_name: fsobral/pack-u


text:
C    @(#)kmpagc.f	20.3 2/13/96
        function kmpagc (m,n)
C
      include 'ipfinc/parametr.inc'

      include 'ipfinc/agc.inc'
      include 'ipfinc/qsdup.inc'
C
        kmpagc = kagc(1,m) - kagc(1,n)
        if (kmpagc .eq. 0 .and. m .ne. n) dupsw = .true.
        return
        end