! { dg-do run }
!
! Check fix for correctly deep copying allocatable components.
! PR fortran/59678
! Contributed by Andre Vehreschild  <vehre@gmx.de>
!
program alloc_comp_copy_test

  type InnerT
    integer :: ii
    integer, allocatable :: ai
    integer, allocatable :: v(:)
  end type InnerT

  type T
    integer :: i
    integer, allocatable :: a_i
    type(InnerT), allocatable :: it
    type(InnerT), allocatable :: vec(:)
  end type T

  type(T) :: o1, o2
  class(T), allocatable :: o3, o4
  o1%i = 42

  call copyO(o1, o2)
  if (o2%i /= 42) call abort ()
  if (allocated(o2%a_i)) call abort()
  if (allocated(o2%it)) call abort()
  if (allocated(o2%vec)) call abort()

  allocate (o1%a_i, source=2)
  call copyO(o1, o2)
  if (o2%i /= 42) call abort ()
  if (.not. allocated(o2%a_i)) call abort()
  if (o2%a_i /= 2) call abort()
  if (allocated(o2%it)) call abort()
  if (allocated(o2%vec)) call abort()

  allocate (o1%it)
  o1%it%ii = 3
  call copyO(o1, o2)
  if (o2%i /= 42) call abort ()
  if (.not. allocated(o2%a_i)) call abort()
  if (o2%a_i /= 2) call abort()
  if (.not. allocated(o2%it)) call abort()
  if (o2%it%ii /= 3) call abort()
  if (allocated(o2%it%ai)) call abort()
  if (allocated(o2%it%v)) call abort()
  if (allocated(o2%vec)) call abort()

  allocate (o1%it%ai)
  o1%it%ai = 4
  call copyO(o1, o2)
  if (o2%i /= 42) call abort ()
  if (.not. allocated(o2%a_i)) call abort()
  if (o2%a_i /= 2) call abort()
  if (.not. allocated(o2%it)) call abort()
  if (o2%it%ii /= 3) call abort()
  if (.not. allocated(o2%it%ai)) call abort()
  if (o2%it%ai /= 4) call abort()
  if (allocated(o2%it%v)) call abort()
  if (allocated(o2%vec)) call abort()

  allocate (o1%it%v(3), source= 5)
  call copyO(o1, o2)
  if (o2%i /= 42) call abort ()
  if (.not. allocated(o2%a_i)) call abort()
  if (o2%a_i /= 2) call abort()
  if (.not. allocated(o2%it)) call abort()
  if (o2%it%ii /= 3) call abort()
  if (.not. allocated(o2%it%ai)) call abort()
  if (o2%it%ai /= 4) call abort()
  if (.not. allocated(o2%it%v)) call abort()
  if (any (o2%it%v /= 5) .or. size (o2%it%v) /= 3) call abort()
  if (allocated(o2%vec)) call abort()

  allocate (o1%vec(2))
  o1%vec(:)%ii = 6
  call copyO(o1, o2)
  if (o2%i /= 42) call abort ()
  if (.not. allocated(o2%a_i)) call abort()
  if (o2%a_i /= 2) call abort()
  if (.not. allocated(o2%it)) call abort()
  if (o2%it%ii /= 3) call abort()
  if (.not. allocated(o2%it%ai)) call abort()
  if (o2%it%ai /= 4) call abort()
  if (.not. allocated(o2%it%v)) call abort()
  if (size (o2%it%v) /= 3) call abort()
  if (any (o2%it%v /= 5)) call abort()
  if (.not. allocated(o2%vec)) call abort()
  if (size(o2%vec) /= 2) call abort()
  if (any(o2%vec(:)%ii /= 6)) call abort()
  if (allocated(o2%vec(1)%ai) .or. allocated(o2%vec(2)%ai)) call abort()
  if (allocated(o2%vec(1)%v) .or. allocated(o2%vec(2)%v)) call abort()

  allocate (o1%vec(2)%ai)
  o1%vec(2)%ai = 7
  call copyO(o1, o2)
  if (o2%i /= 42) call abort ()
  if (.not. allocated(o2%a_i)) call abort()
  if (o2%a_i /= 2) call abort()
  if (.not. allocated(o2%it)) call abort()
  if (o2%it%ii /= 3) call abort()
  if (.not. allocated(o2%it%ai)) call abort()
  if (o2%it%ai /= 4) call abort()
  if (.not. allocated(o2%it%v)) call abort()
  if (size (o2%it%v) /= 3) call abort()
  if (any (o2%it%v /= 5)) call abort()
  if (.not. allocated(o2%vec)) call abort()
  if (size(o2%vec) /= 2) call abort()
  if (any(o2%vec(:)%ii /= 6)) call abort()
  if (allocated(o2%vec(1)%ai)) call abort()
  if (.not. allocated(o2%vec(2)%ai)) call abort()
  if (o2%vec(2)%ai /= 7) call abort()
  if (allocated(o2%vec(1)%v) .or. allocated(o2%vec(2)%v)) call abort()

  allocate (o1%vec(1)%v(3))
  o1%vec(1)%v = [8, 9, 10]
  call copyO(o1, o2)
  if (o2%i /= 42) call abort ()
  if (.not. allocated(o2%a_i)) call abort()
  if (o2%a_i /= 2) call abort()
  if (.not. allocated(o2%it)) call abort()
  if (o2%it%ii /= 3) call abort()
  if (.not. allocated(o2%it%ai)) call abort()
  if (o2%it%ai /= 4) call abort()
  if (.not. allocated(o2%it%v)) call abort()
  if (size (o2%it%v) /= 3) call abort()
  if (any (o2%it%v /= 5)) call abort()
  if (.not. allocated(o2%vec)) call abort()
  if (size(o2%vec) /= 2) call abort()
  if (any(o2%vec(:)%ii /= 6)) call abort()
  if (allocated(o2%vec(1)%ai)) call abort()
  if (.not. allocated(o2%vec(2)%ai)) call abort()
  if (o2%vec(2)%ai /= 7) call abort()
  if (.not. allocated(o2%vec(1)%v)) call abort()
  if (any (o2%vec(1)%v /= [8,9,10])) call abort()
  if (allocated(o2%vec(2)%v)) call abort()

  ! Now all the above for class objects.
  allocate (o3, o4)
  o3%i = 42

  call copyO(o3, o4)
  if (o4%i /= 42) call abort ()
  if (allocated(o4%a_i)) call abort()
  if (allocated(o4%it)) call abort()
  if (allocated(o4%vec)) call abort()

  allocate (o3%a_i, source=2)
  call copyO(o3, o4)
  if (o4%i /= 42) call abort ()
  if (.not. allocated(o4%a_i)) call abort()
  if (o4%a_i /= 2) call abort()
  if (allocated(o4%it)) call abort()
  if (allocated(o4%vec)) call abort()

  allocate (o3%it)
  o3%it%ii = 3
  call copyO(o3, o4)
  if (o4%i /= 42) call abort ()
  if (.not. allocated(o4%a_i)) call abort()
  if (o4%a_i /= 2) call abort()
  if (.not. allocated(o4%it)) call abort()
  if (o4%it%ii /= 3) call abort()
  if (allocated(o4%it%ai)) call abort()
  if (allocated(o4%it%v)) call abort()
  if (allocated(o4%vec)) call abort()

  allocate (o3%it%ai)
  o3%it%ai = 4
  call copyO(o3, o4)
  if (o4%i /= 42) call abort ()
  if (.not. allocated(o4%a_i)) call abort()
  if (o4%a_i /= 2) call abort()
  if (.not. allocated(o4%it)) call abort()
  if (o4%it%ii /= 3) call abort()
  if (.not. allocated(o4%it%ai)) call abort()
  if (o4%it%ai /= 4) call abort()
  if (allocated(o4%it%v)) call abort()
  if (allocated(o4%vec)) call abort()

  allocate (o3%it%v(3), source= 5)
  call copyO(o3, o4)
  if (o4%i /= 42) call abort ()
  if (.not. allocated(o4%a_i)) call abort()
  if (o4%a_i /= 2) call abort()
  if (.not. allocated(o4%it)) call abort()
  if (o4%it%ii /= 3) call abort()
  if (.not. allocated(o4%it%ai)) call abort()
  if (o4%it%ai /= 4) call abort()
  if (.not. allocated(o4%it%v)) call abort()
  if (any (o4%it%v /= 5) .or. size (o4%it%v) /= 3) call abort()
  if (allocated(o4%vec)) call abort()

  allocate (o3%vec(2))
  o3%vec(:)%ii = 6
  call copyO(o3, o4)
  if (o4%i /= 42) call abort ()
  if (.not. allocated(o4%a_i)) call abort()
  if (o4%a_i /= 2) call abort()
  if (.not. allocated(o4%it)) call abort()
  if (o4%it%ii /= 3) call abort()
  if (.not. allocated(o4%it%ai)) call abort()
  if (o4%it%ai /= 4) call abort()
  if (.not. allocated(o4%it%v)) call abort()
  if (size (o4%it%v) /= 3) call abort()
  if (any (o4%it%v /= 5)) call abort()
  if (.not. allocated(o4%vec)) call abort()
  if (size(o4%vec) /= 2) call abort()
  if (any(o4%vec(:)%ii /= 6)) call abort()
  if (allocated(o4%vec(1)%ai) .or. allocated(o4%vec(2)%ai)) call abort()
  if (allocated(o4%vec(1)%v) .or. allocated(o4%vec(2)%v)) call abort()

  allocate (o3%vec(2)%ai)
  o3%vec(2)%ai = 7
  call copyO(o3, o4)
  if (o4%i /= 42) call abort ()
  if (.not. allocated(o4%a_i)) call abort()
  if (o4%a_i /= 2) call abort()
  if (.not. allocated(o4%it)) call abort()
  if (o4%it%ii /= 3) call abort()
  if (.not. allocated(o4%it%ai)) call abort()
  if (o4%it%ai /= 4) call abort()
  if (.not. allocated(o4%it%v)) call abort()
  if (size (o4%it%v) /= 3) call abort()
  if (any (o4%it%v /= 5)) call abort()
  if (.not. allocated(o4%vec)) call abort()
  if (size(o4%vec) /= 2) call abort()
  if (any(o4%vec(:)%ii /= 6)) call abort()
  if (allocated(o4%vec(1)%ai)) call abort()
  if (.not. allocated(o4%vec(2)%ai)) call abort()
  if (o4%vec(2)%ai /= 7) call abort()
  if (allocated(o4%vec(1)%v) .or. allocated(o4%vec(2)%v)) call abort()

  allocate (o3%vec(1)%v(3))
  o3%vec(1)%v = [8, 9, 10]
  call copyO(o3, o4)
  if (o4%i /= 42) call abort ()
  if (.not. allocated(o4%a_i)) call abort()
  if (o4%a_i /= 2) call abort()
  if (.not. allocated(o4%it)) call abort()
  if (o4%it%ii /= 3) call abort()
  if (.not. allocated(o4%it%ai)) call abort()
  if (o4%it%ai /= 4) call abort()
  if (.not. allocated(o4%it%v)) call abort()
  if (size (o4%it%v) /= 3) call abort()
  if (any (o4%it%v /= 5)) call abort()
  if (.not. allocated(o4%vec)) call abort()
  if (size(o4%vec) /= 2) call abort()
  if (any(o4%vec(:)%ii /= 6)) call abort()
  if (allocated(o4%vec(1)%ai)) call abort()
  if (.not. allocated(o4%vec(2)%ai)) call abort()
  if (o4%vec(2)%ai /= 7) call abort()
  if (.not. allocated(o4%vec(1)%v)) call abort()
  if (any (o4%vec(1)%v /= [8,9,10])) call abort()
  if (allocated(o4%vec(2)%v)) call abort()

contains

  subroutine copyO(src, dst)
    type(T), intent(in) :: src
    type(T), intent(out) :: dst

    dst = src
  end subroutine copyO

end program alloc_comp_copy_test


USE DFPORT

character*4 dsaa/'DSAA'/
character*20 scale_line,scale_vpvs,scale,scale_vp,scale_vs
character*2 lv, ver
character*8 ar,md,line
character*1 itt,ps
allocatable dvan(:,:),vvv(:,:),vtmp(:,:)
real hlev(20)
real fia0(100),teta0(100),fib0(100),tetb0(100)
real fmark(200,50),tmark(200,50),smark(200,50)
integer kdot_rgb(3)
integer nmark(200)

common/pi/pi,per
common/center/fi0,tet0
common/keys/key_ft1_xy2

one=1.e0
pi=asin(one)*2.e0
per=pi/180.e0
rz=6371.

open(1,file='../../../model.dat')
read(1,'(a8)')ar
read(1,'(a8)')md
read(1,*)ips
close(1)

write(ps,'(i1)')ips

write(*,*)' ar=',ar,' md=',md,' ps=',ps

i=system('mkdir ..\..\..\TMP_files\att')

!******************************************************************
key_ft1_xy2=1
open(1,file='../../../DATA/'//ar//'/'//md//'/MAJOR_PARAM.DAT')
do i=1,10000
    read(1,'(a8)',end=513)line
    if(line.eq.'GENERAL ') goto 514
end do
513 continue
write(*,*)' cannot find GENERAL INFORMATION in MAJOR_PARAM.DAT!!!'
pause
514 continue
read(1,*)
read(1,*)
read(1,*)
read(1,*)
read(1,*,end=441,err=441)key_ft1_xy2
read(1,*,end=441,err=441)key_true1      ! 1: true locations exist in separate line
read(1,*,end=441,err=441)key_flat1      ! 1: calculations in flat model, 2: spherical velocity
441 close(1)

!******************************************************************
open(1,file='../../../DATA/'//ar//'/'//md//'/MAJOR_PARAM.DAT')
do i=1,10000
	read(1,'(a8)',end=523)line
	if(line.eq.'ATTENUAT') goto 524
end do
523 continue
write(*,*)' cannot find ATTENUATION PARAMETERS in MAJOR_PARAM.DAT!!!'
pause
524 read(1,*)
read(1,*) iter
read(1,*)val11,val12
read(1,*)val21,val22
att_aver=val11; add_att=val12
close(1)

if(ips.eq.2) att_aver=val21; add_att=val22

i=system('mkdir ..\..\..\TMP_files\vert')

key_preview=0
open(1,file='../../../preview_key.txt')
read(1,*,end=771)key_preview
771 close(1)

if(key_preview.ne.0) then
	i=system('mkdir ..\..\..\PICS\'//ar//'\'//md//'\ATT')
557	open(1,file='../../../DATA/'//ar//'/config.txt')
	read(1,*)
	read(1,*) npix_h_x0,npix_h_y0
	read(1,*)
	read(1,*)
	read(1,*) npix_x0,npix_y
	read(1,*)tick_x,tick_y
	read(1,*)
	read(1,*)
	read(1,*)
	read(1,*)
	read(1,*)
	read(1,*)
	read(1,*)scale_line
	read(1,*)a1,a2,dv_min0,dv_max0
	read(1,*)scale_vpvs
	read(1,*)vpvs_min,vpvs_max
	read(1,*,end=559,err=559)scale_vp
	read(1,*,end=559,err=559)vp_min,vp_max
	read(1,*,end=559,err=559)scale_vs
	read(1,*,end=559,err=559)vs_min,vs_max
	close(1)
    goto 558
559 write(*,*)'  '
    write(*,*)' Please update the file "config.txt" in the AREA folder !!!'
    write(*,*)' To the end of this file add lines with color scales for'
    write(*,*)' absolute P and S velocities. For example:'
    write(*,*)
    write(*,*)' rainbow_small.scl	scale for absolute Vp'
    write(*,*)' 2.6 3.2			diapason for absolute Vp'
    write(*,*)' rainbow_small.scl	scale for absolute Vs'
    write(*,*)' 1.7 2.3			diapason for absolute Vs'
    write(*,*)
    write(*,*)' When ready, press ENTER, and the program will make another attempt'
    close(1)
    pause
    goto 557
558 continue

	i=system('copy ..\..\..\COMMON\visual_exe\visual.exe layers.exe')
	i=system('copy ..\..\..\COMMON\scales_scl\'//scale_line//' '//scale_line)
	i=system('copy ..\..\..\COMMON\scales_scl\'//scale_vpvs//' '//scale_vpvs)
	i=system('copy ..\..\..\COMMON\scales_scl\'//scale_vp//' '//scale_vp)
	i=system('copy ..\..\..\COMMON\scales_scl\'//scale_vs//' '//scale_vs)
end if

!******************************************************************
open(1,file='../../../DATA/'//ar//'/'//md//'/MAJOR_PARAM.DAT')
do i=1,10000
	read(1,'(a8)',end=573)line
	if(line.eq.'ORIENTAT') goto 574
end do
573 continue
write(*,*)' cannot find ORIENTATIONS in MAJOR_PARAM.DAT!!!'
pause
574 read(1,*)nornt
close(1)
!******************************************************************

ngr1=1
ngr2=nornt

!******************************************************************
if(key_ft1_xy2.eq.1) then
    open(1,file='../../../DATA/'//ar//'/'//md//'/MAJOR_PARAM.DAT')
    do i=1,10000
        read(1,'(a8)',end=553)line
        if(line.eq.'AREA_CEN') goto 554
    end do
    553 continue
    write(*,*)' cannot find AREA CENTER in MAJOR_PARAM.DAT!!!'
    pause
    554 read(1,*)fi0,tet0
    close(1)
else
    fi0=0
    tet0=0
end if
!******************************************************************

open(1,file='../../../DATA/'//ar//'/'//md//'/MAJOR_PARAM.DAT')
do i=1,10000
	read(1,'(a8)',end=543)line
	if(line.eq.'ATTENUAT') goto 544
end do
543 continue
write(*,*)' cannot find ATTENUATION PARAMETERS in MAJOR_PARAM.DAT!!!'
pause
544 read(1,*)
read(1,*)iter
close(1)

write(itt,'(i1)')iter

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

open(2,file='../../../DATA/'//ar//'/setver.dat')
read(2,*)nver
do ii=1,nver
	read(2,*) fia0(ii),teta0(ii),fib0(ii),tetb0(ii)
end do
read(2,*) dist_from_sec_event
read(2,*) dxsec
read(2,*) zmin,zmax,dzsec
read(2,*) dsmark
read(2,*) smaxx
read(2,*) ismth
551 close(2)

rsmth=ismth+0.5

call read_topo(ar)

do iver=1,nver
    write(ver,'(i2)')iver

    if(key_ft1_xy2.eq.1) then
        fia=fia0(iver)
        teta=teta0(iver)
        fib=fib0(iver)
        tetb=tetb0(iver)
        call SFDEC(fia,teta,0.,xa,ya,Z,fi0,tet0)
        call SFDEC(fib,tetb,0.,xb,yb,Z,fi0,tet0)
    else
        xa=fia0(iver)
        ya=teta0(iver)
        xb=fib0(iver)
        yb=tetb0(iver)
    end if
    write(*,*)' xa=',xa,' ya=',ya
    write(*,*)' xb=',xb,' yb=',yb
    dist=sqrt((xb-xa)*(xb-xa)+(yb-ya)*(yb-ya))
    !write(*,*)' dist=',dist

    sinpov=(yb-ya)/dist
    cospov=(xb-xa)/dist
    nxsec=dist/dxsec+1
    dxsec=dist/(nxsec-1)
    nzsec=(zmax-zmin)/dzsec+1
    dzsec=(zmax-zmin)/(nzsec-1)
    write(*,*)' dist=',dist,' nxsec=',nxsec,' nzsec=',nzsec

    allocate (dvan(nxsec,nzsec),vvv(nxsec,nzsec),vtmp(nxsec,nzsec))
    vvv=0
    dvan=0

    open(11,file='../../../TMP_files/vert/mark_'//ver//'.dat')
    imark=0
    do sss=0.,dist,dsmark
        x=xa+cospov*sss
        y=ya+sinpov*sss
        if(key_ft1_xy2.eq.1) then
            call decsf(x,y,0.,fi0,tet0,FI,TET,h)
        else
            fi=x
            tet=y
        end if
        write(11,*)fi,tet,sss
        imark=imark+1
        fmark(iver,imark)=fi
        tmark(iver,imark)=tet
        smark(iver,imark)=sss
    end do
    imark=imark+1
    fmark(iver,imark)=fib
    tmark(iver,imark)=tetb
    smark(iver,imark)=dist
    close(11)
    nmark(iver)=imark

    ! Draw the position of the section on the surface (line)
    open(11,file='../../../TMP_files/vert/mark_'//ver//'.bln')
    write(11,*) imark
    do i=1,imark
	    write(11,*)fmark(iver,i),tmark(iver,i)	!,smark(i)
    end do
    close(11)

    !Draw topography on the section
    open(11,file='../../../TMP_files/vert/topo_'//ver//'.bln')
    write(11,*)nxsec
    do ix=1,nxsec
        sss=(ix-1)*dxsec
        xcur=xa+((xb-xa)/dist)*sss
        ycur=ya+((yb-ya)/dist)*sss
        if(key_ft1_xy2.eq.1) then
            call decsf(xcur,ycur,0.,fi0,tet0,fff,ttt,h)
        else
            fff=xcur
            ttt=ycur
        end if
        depth=h_lim(fff,ttt)
        ztopo=flat_sph(key_flat1,xcur,ycur,depth)
        write(11,*)sss,-ztopo
    end do
    close(11)

    ! Read the coordinates of the stations
    open(2,file='../../../DATA/'//ar//'/'//md//'/data/stat_xy.dat')
    open(12,file='../../../TMP_files/vert/stat_'//ver//'.dat')
    i=0
    nst1=0
    3   i=i+1
        read(2,*,end=4)xst,yst,depth
        zst = flat_sph(key_flat1,xst,yst,depth)

        xx1=(xst-xa)*cospov+(yst-ya)*sinpov
        yy1=-(xst-xa)*sinpov+(yst-ya)*cospov

        if(abs(yy1).lt.dist_from_sec_event) then
            nst1=nst1+1
            write(12,*)xx1,-zst
        end if
        goto 3
    4   close(2)
    close(12)

    open(1,file='../../../DATA/'//ar//'/'//md//'/data/srces'//itt//'.dat')
    open(11,file='../../../TMP_files/vert/ztr_'//ver//'.dat')
    nzt1=0
    nzt=0
872	read(1,*,end=871)fzt,tzt,depth
        if(key_ft1_xy2.eq.1) then
            call SFDEC(fzt,tzt,0.,xzt,yzt,Z,fi0,tet0)
        else
            xzt=fzt; yzt=tzt
        end if
        zzt = flat_sph(key_flat1,xzt,yzt,depth)
        nzt=nzt+1
        xx1=(xzt-xa)*cospov+(yzt-ya)*sinpov
        yy1=-(xzt-xa)*sinpov+(yzt-ya)*cospov
        if(abs(yy1).lt.dist_from_sec_event) then
            nzt1=nzt1+1
            write(11,*)xx1,-zzt,yy1
        end if
        goto 872
    871 close(1)
    close(11)
    write(*,*)' nst1=',nst1,' nzt1=',nzt1

    vvv=0
    dvan=0

!write(*,*)' dist=',dist,' nxsec=',nxsec

    do igr=ngr1,ngr2
        call prepare_model_att(ar,md,igr,ips)

        do ix=1,nxsec
            sss=(ix-1)*dxsec
            !write(*,*)' ix=',ix,' sss=',sss,' dist=',dist
            !sss=18
            !if(mod(ix,10).eq.0)write(*,*)' ix=',ix,' sss=',sss
            xcur=xa+((xb-xa)/dist)*sss
            ycur=ya+((yb-ya)/dist)*sss

            if(key_ft1_xy2.eq.1) then
                call decsf(xcur,ycur,0.,fi0,tet0,fff,ttt,h)
            else
                fff=xcur; ttt=ycur
            end if
            !write(*,*)' xcur=',xcur,' ycur=',ycur
            !write(*,*)' fi=',fff,' tet=',ttt
            do iz=1,nzsec
                zcur=zmin+(iz-1)*dzsec
                !zcur=5
                depth = sph_flat(key_flat1,xcur,ycur,zcur)

                call att_1_grid(fff,ttt,zcur,smaxx, att,umn)

                depth = h_lim(fff,ttt)
                zlim_up = flat_sph(key_flat1,xcur,ycur,depth)

                if (zcur.lt.zlim_up) umn=0

                !write(*,*)' zcur=',zcur,' att=',att,' umn=',umn

                dvan(ix,iz) = dvan(ix,iz) + att*umn
                vvv(ix,iz)=vvv(ix,iz)+umn
            end do
        end do
    end do

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

    dv_min=9999; dv_max=-9999
    do iz=nzsec,1,-1
        do ix=1,nxsec
            vanom=-999
            if(vvv(ix,iz).gt.0.0001) then
	            van=dvan(ix,iz)/vvv(ix,iz)
                vanom=van*100./att_aver
                if(vanom.lt.dv_min) dv_min=vanom
                if(vanom.gt.dv_max) dv_max=vanom
                !write(*,*)dv_min,dv_max
            end if
            dvan(ix,iz)=vanom
        end do
        !write(*,*)dv_min,dv_max
    end do
write(*,*)' dv_min=',dv_min,' dv_max=',dv_max,' att_aver=',att_aver

! smoothing:
        vtmp=dvan
        do iz=1,nzsec
	        do ix=1,nxsec
		        dvan(ix,iz)=-999
		        if(vvv(ix,iz).lt.0.01) cycle
		        vanm=0.
		        iv=0
		        do ixx=-ismth,ismth
			        if (ix+ixx.lt.1) cycle
			        if (ix+ixx.gt.nxsec) cycle
			        do izz=-ismth,ismth
				        if (iz+izz.lt.1) cycle
				        if (iz+izz.gt.nzsec) cycle
				        if(vvv(ix+ixx,iz+izz).lt.0.01) cycle
				        rr=ixx*ixx+izz*izz
				        r=sqrt(rr)
				        if(r.gt.rsmth) cycle
				        iv=iv+1
				        vanm=vanm+vtmp(ix+ixx,iz+izz)
			        end do
		        end do
		        dvan(ix,iz)=vanm/iv
	        end do
        end do

    open(14,file='../../../TMP_files/att/att_ver'//ps//ver//'.grd')
    write(14,'(a4)')dsaa
    write(14,*)nxsec,nzsec
    write(14,*)0,dist
    write(14,*)-zmax,-zmin
    write(14,*)dv_min,dv_max
    do iz=nzsec,1,-1
	    write(14,*)(dvan(ix,iz),ix=1,nxsec)
    end do
    close(14)

    deallocate (dvan,vvv,vtmp)

!*********************************************************
!*********************************************************
!*********************************************************
!*********************************************************
!*********************************************************
        if(npix_x0.ne.0) then
	        npix_x=npix_x0
        else
	        npix_x = int(npix_y * dist/(zmax-zmin))
        end if
        write(*,*)' section:',ver,' dist=',dist,' npix_x=',npix_x

! PLOT THE ATTENUATION ANOMALIES
        open(14,file='config.txt')
        write(14,*)npix_x,npix_y
        write(14,*)'_______ Size of the picture in pixels (nx,ny)'
        write(14,*)0,dist
        write(14,*)'_______ Physical coordinates along X (xmin,xmax)'
        write(14,*)-zmax,-zmin
        write(14,*)'_______ Physical coordinates along Y (ymin,ymax)'
        write(14,*)tick_x,tick_y
        write(14,*)'_______ Spacing of ticks on axes (dx,dy)'
        write(14,51)ar,md,ps,ver
        51 format('..\..\..\PICS\',a8,'\',a8,'\ATT\ver_att',a1,a2,'.png')
        write(14,*)'_______ Path of the output picture'
        if(ips.eq.1) then
	        write(14,3431)	iver, iver
        3431 format(' P attenuation, section=',i1,'A - ',i1,'B')
        else
	        write(14,3432)	iver,iver
        3432 format(' S attenuation, section=',i1,'A - ',i1,'B')
        end if
        write(14,*)'_______ Title of the plot on the upper axe'
        write(14,*)	1
        write(14,*)'_______ Number of layers'

        59 format('********************************************')

        write(14,59)
        write(14,*)	1
        write(14,*)'_______ Key of the layer (1: contour, 2: line, 3:dots)'
        write(14,52)ps,ver
        52 format('..\..\..\TMP_files\att\att_ver',a1,a2,'.grd')
        write(14,*)'_______ Location of the GRD file'
        write(14,53)scale_line
        53 format(a20)
        write(14,*)'_______ Scale for visualization'
        write(14,*)	dv_min0,dv_max0
        write(14,*)'_______ scale diapason'

        write(14,59)
        write(14,*)	2
        write(14,*)'_______ Key of the layer (1: contour, 2: line, 3:dots)'
        write(14,2143)ver
        2143 format('..\..\..\TMP_files\vert\topo_',a2,'.bln')
        write(14,*)'_______ Location of the BLN file'
        write(14,*)	3
        write(14,*)'_______ Thickness of line in pixels'
        write(14,*)	0,0,0
        write(14,*)'_______ RGB color'

        close(14)

        i=system('layers.exe')

end do

stop
end


text:
! { dg-do run }
! PR35627 Namelist read problem with short logical followed by read real
program test
  implicit none
    LOGICAL :: nlco(200)  ! (1:nbeam)
    REAL(kind=8):: xlbtna(200)  ! (1:nbeam)
  NAMELIST/nbdrive_naml/ nlco, xlbtna
    INTEGER :: nbshapa(200)  ! (1:nbeam)
  NAMELIST/nbdrive_naml/ nbshapa
  nlco = .false.
  xlbtna = 0.0_8
  nbshapa = 0
  open(10, file='t.nml')
  write(10,'(a)') "&nbdrive_naml"
  write(10,'(a)') "nlco = 4*T,"
  write(10,'(a)') "xlbtna = 802.8, 802.8, 802.8, 802.8"
  write(10,'(a)') "nbshapa = 4*1"
  write(10,'(a)') "/"
  rewind(10)
  read(10, nbdrive_naml)
  !write(*,nbdrive_naml)
  close(10, status="delete")
end program test

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
cdis
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis
cdis
cdis
cdis
        subroutine insert_pireps(i4time,cld_hts                          ! I
     1        ,default_clear_cover                                       ! I
     1        ,cld_snd,wt_snd,i_snd,j_snd,n_cld_snd,max_cld_snd          ! I/O
     1        ,lat,lon,ni,nj,nk,ix_low,ix_high,iy_low,iy_high,max_pireps ! I
     1        ,istatus)                                                  ! O

        real lat(ni,nj),lon(ni,nj)

        real cld_hts(nk)

!       Arrays for cloud soundings
        real cld_snd(max_cld_snd,nk)
        real wt_snd(max_cld_snd,nk)
        integer i_snd(max_cld_snd)
        integer j_snd(max_cld_snd)

        character*9  asc9_tim_pirep,asc9_tim_rcvd
        character*80 string

        character*31    ext

        logical l_good_pirep

!       1 Feb 1996     Steve Albers        Remove spread( calls

        max_layers = 3  ! Defined in PIN file and NetCDF RAW data

        lun = 11
        ext = 'pin'
        call open_lapsprd_file_read(lun,i4time,ext,istatus)
        if(istatus .ne. 1)go to 998

        n_cloud = 0
        num_pirep = 0
        num_good_pirep = 0

5       read(11,101,end=900,err=50)string(1:6)
101     format(a6)

50      if(string(2:5) .eq. 'Time')then
!           a9time = string(30:39)
            read(11,151)asc9_tim_pirep,asc9_tim_rcvd
151         format(1x,a9,2x,a9)
            write(6,*)
            write(6,151)asc9_tim_pirep,asc9_tim_rcvd
        endif

        if(string(2:4) .eq. 'Lat')then
            read(11,201)rlat,rlon,ralt
201         format(2(f8.3,2x), f6.0,2i5)
            call latlon_to_rlapsgrid(rlat,rlon,lat,lon,ni,nj,ri,rj
     1                              ,istatus)
!           if(istatus.ne.1)return
            ilaps = nint(ri)
            jlaps = nint(rj)
            write(6,201)rlat,rlon,ralt,ilaps,jlaps
        endif

        if(string(2:5) .eq. 'Clou')then

            if(                   ilaps .ge. ix_low   ! In bounds
     1                      .and. ilaps .le. ix_high  ! In bounds
     1                      .and. jlaps .ge. iy_low   ! In bounds
     1                      .and. jlaps .le. iy_high  ! In bounds
     1                                          )then ! In bounds

              n_cld_snd = n_cld_snd + 1  ! Create a new souding
              num_pirep = num_pirep + 1

              if(num_pirep .gt. max_pireps)then
                  write(6,*)' insert_pireps: Error, too many pireps,'
     1                     ,' check N_PIREP parameter'
                  istatus = 0
                  return
              endif

              l_good_pirep = .false.

!             PIN file input heights are feet MSL
              do i = 1,max_layers
                read(11,203,err=500)cbase_ft,ctop_ft,icover
203             format (12x,2f8.0,i5)
                write(6,*,err=500)' Got a cloud report',
     1                  nint(cbase_ft),nint(ctop_ft),icover

                if(icover .gt. 0)then ! Good Report (icover is in eighths)

                  if(cbase_ft .ge. 0. .and. ctop_ft .ge. 0.)then ! Fill Layer

                    cbase_m = cbase_ft / 3.281
                    ctop_m  = ctop_ft  / 3.281

                    cbuf_low = cbase_m - 500.
                    cbuf_high = ctop_m + 500.

                    cover_rpt = float(icover) / 8.0 ! 1.0

                    n_cloud = n_cloud + 1
                    write(6,*)' Good layer report:      '
     1                       ,nint(cbase_m),nint(ctop_m),cover_rpt

                    do k=1,nk
                        cover = cover_rpt

!                       Fill in Cloud Layer
                        if(cld_hts(k) .ge. cbase_m .and.
     1                     cld_hts(k) .lt. ctop_m                  )then
                            call spread2(cld_snd,wt_snd,i_snd,j_snd
     1                                  ,n_cld_snd,max_cld_snd
     1                                  ,nk,ilaps,jlaps,k,cover,1.)
                            write(6,*)' Fill in k = ',k,cover,cld_hts(k)
                            l_good_pirep = .true.
                        endif

                        cover = default_clear_cover

!                       Fill in clear buffer under cloud layer
                        if(cld_hts(k) .ge. cbuf_low .and.
     1                     cld_hts(k) .lt. cbase_m                 )then
                            call spread2(cld_snd,wt_snd,i_snd,j_snd
     1                                  ,n_cld_snd,max_cld_snd
     1                                  ,nk,ilaps,jlaps,k,cover,1.)
                            write(6,*)' Fill in k = ',k,cover,cld_hts(k)
                            l_good_pirep = .true.
                        endif

!                       Fill in clear buffer above cloud layer
                        if(cld_hts(k) .ge. ctop_m .and.
     1                     cld_hts(k) .lt. cbuf_high               )then
                            call spread2(cld_snd,wt_snd,i_snd,j_snd
     1                                  ,n_cld_snd,max_cld_snd
     1                                  ,nk,ilaps,jlaps,k,cover,1.)
                            write(6,*)' Fill in k = ',k,cover,cld_hts(k)
                            l_good_pirep = .true.
                        endif

                    enddo ! K cld_hts

                  elseif(cbase_ft .gt. 0)then
                    cbase_m = cbase_ft / 3.281
                    ctop_m  = cbase_m + 1000.

                    cbuf_low = cbase_m - 500.

                    cover_rpt = float(icover) / 8.0 ! 1.0

                    write(6,*)' Only a base reported:'
     1                       ,nint(cbase_m),nint(ctop_m),cover_rpt

                    n_cloud = n_cloud + 1

                    do k=1,nk
                        cover = cover_rpt

!                       Fill in Cloud Layer
                        if(cld_hts(k) .ge. cbase_m .and.
     1                     cld_hts(k) .lt. ctop_m                  )then
                            call spread2(cld_snd,wt_snd,i_snd,j_snd
     1                                  ,n_cld_snd,max_cld_snd
     1                                  ,nk,ilaps,jlaps,k,cover,1.)
                            write(6,*)' Fill in k = ',k,cover,cld_hts(k)
                            l_good_pirep = .true.
                        endif

                        cover = default_clear_cover

!                       Fill in clear buffer under cloud layer
                        if(cld_hts(k) .ge. cbuf_low .and.
     1                     cld_hts(k) .lt. cbase_m                 )then
                            call spread2(cld_snd,wt_snd,i_snd,j_snd
     1                                  ,n_cld_snd,max_cld_snd
     1                                  ,nk,ilaps,jlaps,k,cover,1.)
                            write(6,*)' Fill in k = ',k,cover,cld_hts(k)
                            l_good_pirep = .true.
                        endif

                    enddo ! K cld_hts

                  elseif(ctop_ft .gt. 0)then
                    ctop_m  = ctop_ft  / 3.281
                    cbase_m = ctop_m - 500.

                    cbuf_high = ctop_m + 500.

                    cover_rpt = float(icover) / 8.0 ! 1.0

                    write(6,*)' Only a top reported:'
     1                       ,nint(cbase_m),nint(ctop_m),cover_rpt

                    n_cloud = n_cloud + 1

                    do k=1,nk
                        cover = cover_rpt

!                       Fill in Cloud Layer
                        if(cld_hts(k) .ge. cbase_m .and.
     1                     cld_hts(k) .lt. ctop_m                  )then
                            call spread2(cld_snd,wt_snd,i_snd,j_snd
     1                                  ,n_cld_snd,max_cld_snd
     1                                  ,nk,ilaps,jlaps,k,cover,1.)
                            write(6,*)' Fill in k = ',k,cover,cld_hts(k)
                            l_good_pirep = .true.
                        endif

                        cover = default_clear_cover

!                       Fill in clear buffer above cloud layer
                        if(cld_hts(k) .ge. ctop_m .and.
     1                     cld_hts(k) .lt. cbuf_high               )then
                            call spread2(cld_snd,wt_snd,i_snd,j_snd
     1                                  ,n_cld_snd,max_cld_snd
     1                                  ,nk,ilaps,jlaps,k,cover,1.)
                            write(6,*)' Fill in k = ',k,cover,cld_hts(k)
                            l_good_pirep = .true.
                        endif

                    enddo ! K cld_hts

                  endif ! Base and/or Top reported

                endif ! Cover > 0

              enddo ! i cloud layer

              if(l_good_pirep)num_good_pirep = num_good_pirep + 1

            else ! Out of Bounds
!             Skip over cloud report
              do i = 1,max_layers
                read(11,203,err=500)cbase_ft,ctop_ft,icover
              enddo

              write(6,*)' Out of domain perimeter'

            endif ! In Bounds

        endif ! Cloud Report String

!       if(string(2:5) .eq. 'Sky ')then
!           read(11,204,err=500)isky_cover
!204         format (40x,i4)
!            write(6,*)' sky cover = ',isky_cover
!        endif

500     goto5

900     write(6,*)
        write(6,*)' Completed insertion of pireps'
        write(6,*)' num pireps/num good pireps/cloud layers = '
     1                  ,num_pirep,num_good_pirep,n_cloud
        close(11)

        istatus = 1
        return

998     write(6,*)' Warning, could not find the pirep file'
        istatus = 1
        return

        end

ext: f

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


! { dg-do compile }
! Tests the fix for PR55984.
!
! Contributed by Sylwester Arabas  <slayoo@staszic.waw.pl>
!
module bcd_m
  type, abstract :: bcd_t
    contains
    procedure(bcd_fill_halos), deferred :: fill_halos
  end type
  abstract interface
    subroutine bcd_fill_halos(this)
      import :: bcd_t
      class(bcd_t ) :: this
    end subroutine
  end interface
end module

module solver_m
  use bcd_m
  type, abstract :: solver_t
    integer :: n, hlo
    class(bcd_t), pointer :: bcx, bcy
    contains
    procedure(solver_advop), deferred :: advop
  end type
  abstract interface
    subroutine solver_advop(this)
      import solver_t
      class(solver_t) :: this
    end subroutine
  end interface
  contains
end module

module solver_mpdata_m
  use solver_m
  type :: mpdata_t
    class(bcd_t), pointer :: bcx, bcy
    contains
    procedure :: advop => mpdata_advop
  end type
  contains
  subroutine mpdata_advop(this)
    class(mpdata_t) :: this
    associate ( bcx => this%bcx, bcy => this%bcy )
      call bcx%fill_halos()
    end associate
  end subroutine
end module

  use solver_mpdata_m
  class(mpdata_t), allocatable :: that
  call mpdata_advop (that)
end


text:
\ primitive circular buffer SwiftForth                             uh 2021-08-05

\ see https://github.com/Project-Forth-Works/data-structures/primitive-circular-buffer for details

\ Internal words
package primitive-circular-buffer-package

\ 1. define the fields of the primitive circular buffer { k idx x0 ... xk-1 }
0
      : _idx ( -- ) @ ;                  \ index to next free entry
cell+ : _x0 ( addr1 -- addr2 ) @ cell+ ; \ address of first entry
1 cells +Field _k                        \ number of entries in primitive circular buffer
drop

\ 2. access entries in primitive circular buffer pcb
: _x[] ( i pcb -- 'x[i]) \ get address of ith entry in primitive circular buffer pcb
    _x0  swap cells + ;

: _x[idx]  ( pcb -- 'x[idx] ) \ get address of next free entry in primitive circular pcb
    dup  _idx @  swap _x[] ;

\ 3. manage circular buffer index
: advance ( k i -- i' ) \ increment index i wrapping around with k entries
   1+ swap mod ;

: advance-idx ( pcb -- ) \ increment index of primitive circular buffer pcb
   dup _k @  swap _idx dup >r  @ advance r> ! ;

\ 4. Interface words
public

: pcb-init ( pcb -- )  \ initialize primitive circular buffer
   0 swap _idx ! ;

: pcb! ( x pcb -- addr )  \ store value x in next free entry of primitive circular buffer cb
   swap over _x[idx] !  advance-idx ;

: pcb@0 ( pcb -- x ) \ retrieve the oldest entry from primitive circular buffer cb
   _x[idx] @ ;

: pcb@ ( i pcb -- x ) \ retrieve the ith oldest entry from primitive circular buffer cb
   dup >r _idx @ +  r@ _k @ mod  r> _x[] @ ;

private

: allocate-primitive-circular-buffer ( k -- addr ) \ allocate a primitive circular buffer leaving its address
   here swap cells cell+ ( for idx ) allot ; \ RAM: ( idx x0 ... x_k-1 )

: setup-primitive-circular-buffer ( k addr -- )
   , ( 'ram ) , ( k ) ; \ ROM: { 'ram k }

public

: new-primitive-circular-buffer ( k -- addr )
   dup allocate-primitive-circular-buffer \ RAM
   setup-primitive-circular-buffer ;      \ ROM

: Primitive-Circular-Buffer ( k <name> -- )
   dup allocate-primitive-circular-buffer \ RAM
   Create
   setup-primitive-circular-buffer ; \ ROM

end-package

\ 5. integrated test

Marker ***tests***

include %SwiftForth/unsupported/anstest/tester.f

5 Primitive-Circular-Buffer pcb  pcb pcb-init

{> 10 pcb pcb! 20 pcb pcb! 30 pcb pcb!
   40 pcb pcb! 50 pcb pcb! 60 pcb pcb!
    0 pcb pcb@ 70 pcb pcb!  0 pcb pcb@  0 pcb pcb!  0 pcb pcb@  0 pcb pcb!
    0 pcb pcb@  0 pcb pcb!  0 pcb pcb@  0 pcb pcb!  0 pcb pcb@  0 pcb pcb!
    -> 20 30 40 50 60 70 }

***tests***

max_stars_repo_name: embeddingforth/embeddinForth

max_issues_repo_name: embeddingforth/embeddinForth

max_forks_repo_name: embeddingforth/embeddinForth


! enum_space.f90 --
!     Module with subroutines to enumerate grid points in space
!
!     Motivation:
!     When solving Diophantine equations with brute force (so simply try
!     all combinations of integer values and see whether that gives a
!     solution), one would like to search a limited region of the solution
!     space and stop after, say, 100 solutions. For instance:
!
!     x**2 + y**2 = z**3 + 3*z**2
!
!     You never know in which region these solutions lie, so it may that
!     you need to examine a very large range for all variables involved:
!
!     do z = 1,1000
!         do y = 1,1000
!             do x = 1,1000
!                 ... is (x,y,z) a solution? Print it
!             enddo
!         enddo
!     enddo
!
!     The drawback is that you examine the cube line by line and if
!     the solutions have coordinates of the same order of magnitude,
!     So examining solutions with very different sizes (x=1000, y= 10, z=1)
!     may waste a lot of time.
!     The procedures in this module use a simple enumeration technique that
!     makes it possible to enumerate the grid points around the origin
!     in ever wider shells (square around the origin or an octaeder in 3D
!     space):
!
!     |
!     9
!     |
!     5   8   .
!     | .   .   .
!     2   4   7   .
!     | .   .   .   .
!     0---1---3---6--10-
!
!     This works for all dimensions - just repeatedly split up the
!     coordinates (see the implementation of enum_octant), but it
!     is non-trivial to expand this to an arbitrary dimension.
!
!     To cover both positive and negative integers, the module uses an
!     alternating scheme: 0, 1, -1, 2, -2, 3, -3, ...
!
!     TODO: enum_ndspace? Non-trivial as said
!
module enumerate_space
    implicit none

    integer, parameter :: start_enumeration = 0

contains

! enum_1dspace --
!     Enumerate the grid points on a line
!
! Arguments:
!     index         The index of the grid point (updated to give the next
!                   index, unless you specify update = .false.)
!     x             The x-coordinate of the grid point (output)
!     update        (Optional) whether to update the index (default) or not
!
! Result:
!     The x-coordinate of the grid point with the given index
!
! Note:
!     Initialise the index to "start_enumeration"
!
subroutine enum_1dspace( index, x, update )
    integer, intent(inout)        :: index
    integer, intent(out)          :: x
    logical, intent(in), optional :: update

    logical                       :: do_update

    do_update = .true.
    if ( present(update) ) then
        do_update = update
    endif

    x = (index + 1) / 2
    if ( mod(index,2) == 0 ) then
        x = -x
    endif

    if ( do_update ) then
        index = index + 1
    endif
end subroutine enum_1dspace

! enum_2dspace --
!     Enumerate the grid points on a plane
!
! Arguments:
!     index         The index of the grid point (updated to give the next
!                   index, unless you specify update = .false.)
!     x             The x-coordinate of the grid point (output)
!     y             The y-coordinate of the grid point (output)
!     update        (Optional) whether to update the index (default) or not
!
! Result:
!     The x- and y-coordinates of the grid point with the given index
!
! Note:
!     Initialise the index to "start_enumeration"
!
subroutine enum_2dspace( index, x, y, update )
    integer, intent(inout)        :: index
    integer, intent(out)          :: x
    integer, intent(out)          :: y
    logical, intent(in), optional :: update

    logical                       :: do_update
    integer                       :: xx, yy

    do_update = .true.
    if ( present(update) ) then
        do_update = update
    endif

    call enum_quadrant( index, xx, yy, do_update )
    call enum_1dspace( xx, x )
    call enum_1dspace( yy, y )

end subroutine enum_2dspace

! enum_3dspace --
!     Enumerate the grid points in 3D space
!
! Arguments:
!     index         The index of the grid point (updated to give the next
!                   index, unless you specify update = .false.)
!     x             The x-coordinate of the grid point (output)
!     y             The y-coordinate of the grid point (output)
!     z             The z-coordinate of the grid point (output)
!     update        (Optional) whether to update the index (default) or not
!
! Result:
!     The x- and y-coordinates of the grid point with the given index
!
! Note:
!     Initialise the index to "start_enumeration"
!
subroutine enum_3dspace( index, x, y, z, update )
    integer, intent(inout)        :: index
    integer, intent(out)          :: x
    integer, intent(out)          :: y
    integer, intent(out)          :: z
    logical, intent(in), optional :: update

    logical                       :: do_update
    integer                       :: xx, yy, zz

    do_update = .true.
    if ( present(update) ) then
        do_update = update
    endif

    call enum_octant( index, xx, yy, zz, do_update )
    call enum_1dspace( xx, x )
    call enum_1dspace( yy, y )
    call enum_1dspace( zz, z )

end subroutine enum_3dspace

! enum_quadrant --
!     Enumerate the grid points in the first 2D quadrant
!
! Arguments:
!     index         The index of the grid point (updated to give the next
!                   index, unless you specify update = .false.)
!     x             The x-coordinate of the grid point (output)
!     y             The y-coordinate of the grid point (output)
!     update        (Optional) whether to update the index (default) or not
!
! Result:
!     The x- and y-coordinates of the grid point with the given index
!
! Note:
!     Initialise the index to "start_enumeration"
!
subroutine enum_quadrant( index, x, y, update )
    integer, intent(inout)        :: index
    integer, intent(out)          :: x
    integer, intent(out)          :: y
    logical, intent(in), optional :: update

    integer, parameter            :: dp = kind(1.0d0)
    real(kind=dp)                 :: dindex
    real(kind=dp)                 :: dn
    integer                       :: n
    logical                       :: do_update

    do_update = .true.
    if ( present(update) ) then
        do_update = update
    endif

    !
    ! Use double-precision reals to avoid rounding and overflow problems
    ! And take care of very large indices
    !
    dindex = index

    dn     = -0.5_dp + sqrt( 200.0 * dindex + 25.0_dp + 1.0_dp ) / 10.0_dp
    n      = dn

    if ( mod(n,2) == 0 ) then
        y = index - (n / 2) * (n + 1)
    else
        y = index - n * ((n + 1) / 2)
    endif
    x  = n - y

    if ( do_update ) then
        index = index + 1
    endif
end subroutine enum_quadrant

! enum_octant --
!     Enumerate the grid points in the first 3D octant
!
! Arguments:
!     index         The index of the grid point (updated to give the next
!                   index, unless you specify update = .false.)
!     x             The x-coordinate of the grid point (output)
!     y             The y-coordinate of the grid point (output)
!     z             The y-coordinate of the grid point (output)
!     update        (Optional) whether to update the index (default) or not
!
! Result:
!     The x-, y and z-coordinates of the grid point with the given index
!
! Note:
!     Initialise the index to "start_enumeration"
!
subroutine enum_octant( index, x, y, z, update )
    integer, intent(inout)        :: index
    integer, intent(out)          :: x
    integer, intent(out)          :: y
    integer, intent(out)          :: z
    logical, intent(in), optional :: update

    integer, parameter            :: dp = kind(1.0d0)
    real(kind=dp)                 :: dindex
    real(kind=dp)                 :: dn
    integer                       :: n
    integer                       :: nindex, prev_index
    logical                       :: do_update

    do_update = .true.
    if ( present(update) ) then
        do_update = update
    endif

    !
    ! We can not use the same approach as for 2D, as it difficult to solve
    ! the third-degree equation that would be involved. Use the brute force
    ! approach instead.
    !
    ! For each plane we have n(n+1)/2 points, so count how many such planes
    ! are completely filled. Some manipulation of the numbers involved is
    ! required to make sure all points are covered.
    !
    nindex = index + 1
    n      = 1

    do while( nindex > 0 )
        prev_index = nindex
        if ( mod(n,2) == 0 ) then
            nindex = nindex - (n / 2) * (n + 1)
        else
            nindex = nindex - n * ((n + 1) / 2)
        endif

        n = n + 1
    enddo

    nindex = prev_index - 1
    n      = n - 2
    call enum_quadrant( nindex, x, y )
    z = n - x - y

    if ( do_update ) then
        index = index + 1
    endif
end subroutine enum_octant
end module enumerate_space


text: Prestige Equestrians offers English Hunter and Jumper Horses horse sales, riding lessons and training in the Davis California area at Hacienda Halguena.  We specialized in Hunter Jumper Training and offer a wide variety of horse riding lesson and training options.  Feel free to leave us a message with any questions and Thanks for looking at our profile!

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
MODULE META
!==================================
! Global variables for metadynamics
!==================================
  IMPLICIT NONE

  ! Input Parameters
  integer :: m_pwr, n_pwr, num_cv
  integer :: m2n_pwr, n2n_pwr, msf_pwr, nsf_pwr
  integer :: metaperiod, NG, metaframe
  real(8) :: gheight, gwidth, g2nwidth, gsfwidth
  real(8) :: rzero               ! parameter of NN distance distribution [bulk latt ref]
  real(8) :: r2n, rsf, d2n, dsf  ! parameters [bulk latt ref]
  logical :: collvar_wanted      ! to calculate CV also during a microcan
  INTEGER :: cn_n_pwr, cn_m_pwr  ! cn bim
  REAL(8) :: cn_rzero, cn_gwidth ! cn bim
  LOGICAL :: cn_aa, cn_bb, cn_ab ! cn bim
  LOGICAL :: yesno_wanted        ! cn bim
  REAL(8) :: d_com_gwidth        ! d_com
  INTEGER :: natoma, natomb      ! d_com
  LOGICAL :: d_com_wanted        ! d_com
  character(LEN=20), dimension(2) :: collvar_name ! deps on number of cv -- max 2
  integer, dimension(2)           :: collvar_case ! to associate a CV to a number (for select statement)

  ! Common Neighbour Number
  REAL(8) :: cnf_num, gcnnwidth, inflim, suplim
  REAL(8), DIMENSION(14):: cnf_in
  REAL(8), DIMENSION(5) :: win_in

  ! arrays used as local, but defined global not to reallocate it
  double precision, allocatable :: deno (:), halfdeno(:)
  double precision, allocatable :: dS(:,:)
  double precision, allocatable :: dVg_ds(:)
  double precision, allocatable :: expoc(:)

  double precision, allocatable :: gausswidth(:)
  double precision, ALLOCATABLE :: s_of_t(:,:)
  double precision, ALLOCATABLE :: coll_var(:)   ! vector of current values of collective variables (num_cv)

  ! Meta-force Computation CN
  double precision, ALLOCATABLE :: s(:,:)        ! s(atom-i, atom-j)
  double precision, ALLOCATABLE :: CNatom(:)     ! coord. no. of every atoms
  double precision, ALLOCATABLE :: fgx(:),fgy(:),fgz(:) ! to natom
  DOUBLE PRECISION, ALLOCATABLE :: dS_dx(:), dS_dy(:), dS_dz(:) ! to natom

  ! Meta-force Computation C2N
  DOUBLE PRECISION, ALLOCATABLE :: s2n(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: C2Natom(:)
  DOUBLE PRECISION, ALLOCATABLE :: fg2nx(:), fg2ny(:), fg2nz(:)
  DOUBLE PRECISION, ALLOCATABLE :: dS2n_dx(:), dS2n_dy(:), dS2n_dz(:)

  ! Meta-force Computation SFN
  DOUBLE PRECISION, ALLOCATABLE :: ssf(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: SFNatom(:)
  DOUBLE PRECISION, ALLOCATABLE :: fgsfx(:), fgsfy(:), fgsfz(:)
  DOUBLE PRECISION, ALLOCATABLE :: dSsf_dx(:), dSsf_dy(:), dSsf_dz(:)

  ! Meta-force Computation CNF
  DOUBLE PRECISION, ALLOCATABLE :: fgcnnx(:), fgcnny(:), fgcnnz(:)
  DOUBLE PRECISION, ALLOCATABLE :: dScnn_dx(:), dScnn_dy(:), dScnn_dz(:)

  ! Meta-force Computation CN_bim
  double precision, ALLOCATABLE :: cn_s(:,:)        ! s(atom-i, atom-j)
  double precision, ALLOCATABLE :: cn_CNatom(:)     ! coord. no. of every atoms
  double precision, ALLOCATABLE :: cn_fgx(:),cn_fgy(:),cn_fgz(:) ! to natom
  DOUBLE PRECISION, ALLOCATABLE :: cn_dS_dx(:), cn_dS_dy(:), cn_dS_dz(:) ! to natom
  LOGICAL, ALLOCATABLE, DIMENSION (:,:) :: pair_yesno_mat

  ! Meta-force Computation d_com
  DOUBLE PRECISION              :: d_com_cv
  DOUBLE PRECISION, ALLOCATABLE :: d_com_dcv_dx(:), d_com_dcv_dy(:), d_com_dcv_dz(:)
  DOUBLE PRECISION, ALLOCATABLE :: d_com_fgx(:), d_com_fgy(:), d_com_fgz(:)

  CONTAINS

  !----------------------------------------------------------------------------
  ! Sigmoid function with and without derivatives
  !----------------------------------------------------------------------------
  ELEMENTAL FUNCTION f_sigmoid(dummy, pairkind)
     ! This function is elemental because is called using array arguments
     ! inside a WHERE statement (in comneighfun_no_light)

     ! cnf_in(1)  = d0
     ! cnf_in(2)  = r0
     ! cnf_in(3)  = c0
     ! cnf_in(4)  = m
     ! cnf_in(5)  = n
     ! cnf_in(6)  = d0 A-A pairkind = 1
     ! cnf_in(7)  = d0 B-B pairkind = 2
     ! cnf_in(8)  = d0 A-B pairkind = 3
     ! cnf_in(9)  = r0 A-A pairkind = 1
     ! cnf_in(10) = r0 B-B pairkind = 2
     ! cnf_in(11) = r0 A-B pairkind = 3
     ! cnf_in(12) = c0 A-A pairkind = 1
     ! cnf_in(13) = c0 B-B pairkind = 2
     ! cnf_in(14) = c0 A-B pairkind = 3

     ! cnf_calc(1) = rij
     ! cnf_calc(2) = rij-d0
     ! cnf_calc(3) = h0
     ! cnf_calc(4) = (rij-d0)/r0
     ! cnf_calc(4) = (rij-d0-c0)/(r0-c0)

     IMPLICIT NONE
     REAL(8) :: f_sigmoid
     REAL(8), DIMENSION(4) :: cnf_calc
     REAL(8), INTENT(IN) :: dummy
     INTEGER(4), INTENT(IN) :: pairkind

     cnf_calc(1) = dummy
     cnf_calc(2) = cnf_calc(1) - cnf_in(pairkind+5) ! cnf_calc(2) = (rij - d0)

     IF (cnf_calc(2).LE.0.d0) THEN
        f_sigmoid = 1.d0
     ELSE
        cnf_calc(3) = 1/(1-(cnf_in(5)*cnf_in(pairkind+8)/ &
                   & ( cnf_in(4)*(cnf_in(pairkind+8)-cnf_in(pairkind+11)) )) )

        IF (    (0.d0.LT.cnf_calc(2)) .AND.(cnf_calc(2).LE.cnf_in(pairkind+8))) THEN
           cnf_calc(4) = cnf_calc(2)/cnf_in(pairkind+8)
           f_sigmoid = (cnf_calc(3)-1.0d0)* cnf_calc(4)**cnf_in(4) + 1.0d0

        ELSE  ! (Note: this will work in combination with a cutoff matrix)
           cnf_calc(4) = (cnf_calc(2)-cnf_in(pairkind+11))/(cnf_in(pairkind+8)-cnf_in(pairkind+11))
           f_sigmoid = cnf_calc(3)* cnf_calc(4)**cnf_in(5)

        END IF
     END IF
  END FUNCTION f_sigmoid

  FUNCTION f_sigmoid_der(dummy, pairkind)

     ! cnf_in(1)  = d0
     ! cnf_in(2)  = r0
     ! cnf_in(3)  = c0
     ! cnf_in(4)  = m
     ! cnf_in(5)  = n
     ! cnf_in(6)  = d0 A-A pairkind = 1
     ! cnf_in(7)  = d0 B-B pairkind = 2
     ! cnf_in(8)  = d0 A-B pairkind = 3
     ! cnf_in(9)  = r0 A-A pairkind = 1
     ! cnf_in(10) = r0 B-B pairkind = 2
     ! cnf_in(11) = r0 A-B pairkind = 3
     ! cnf_in(12) = c0 A-A pairkind = 1
     ! cnf_in(13) = c0 B-B pairkind = 2
     ! cnf_in(14) = c0 A-B pairkind = 3

     ! cnf_calc(1) = rij
     ! cnf_calc(2) = rij-d0
     ! cnf_calc(3) = h0
     ! cnf_calc(4) = (rij-d0)/r0
     ! cnf_calc(4) = (rij-d0-c0)/(r0-c0)
     ! cnf_calc(5) = (h0-1) * ((rij-d0)/r0)**(m-1)
     ! cnf_calc(5) = h0 * ((rij-d0-c0)/(r0-c0))**(n-1)

     ! f_sigmoid_der(1) = f_sigmoid
     ! f_sigmoid_der(2) = derivative of f_sigmoid with respect to the pair distance rij

     IMPLICIT NONE
     REAL(8), DIMENSION(2) :: f_sigmoid_der
     REAL(8), DIMENSION(5) :: cnf_calc
     REAL(8)               :: dummy
     INTEGER(4)            :: pairkind

     cnf_calc(1) = dummy
     cnf_calc(2) = cnf_calc(1) - cnf_in(pairkind+5) ! cnf_calc(2) = (rij - d0)

     IF (cnf_calc(2).LE.0.d0) THEN
        f_sigmoid_der(1) = 1.d0
        f_sigmoid_der(2) = 0.d0
     ELSE
        cnf_calc(3) = 1/(1-(cnf_in(5)*cnf_in(pairkind+8)/ &
                   & ( cnf_in(4)*(cnf_in(pairkind+8)-cnf_in(pairkind+11)) )) )

        IF (    (0.d0.LT.cnf_calc(2)) .AND.(cnf_calc(2).LE.cnf_in(pairkind+8))) THEN
           cnf_calc(4) = cnf_calc(2)/cnf_in(pairkind+8)
           cnf_calc(5) = (cnf_calc(3)-1.d0) * cnf_calc(4)**(cnf_in(4)-1)

           f_sigmoid_der(1) = cnf_calc(4) * cnf_calc(5) +1.0d0
           f_sigmoid_der(2) = cnf_in(4) * (1/cnf_in(pairkind+8)) * cnf_calc(5)

        ELSE   ! (Note: this will work in combination with a cutoff matrix)
           cnf_calc(4) = (cnf_calc(2)-cnf_in(pairkind+11))/(cnf_in(pairkind+8)-cnf_in(pairkind+11))
           cnf_calc(5) = cnf_calc(3) * cnf_calc(4)**(cnf_in(5)-1)

           f_sigmoid_der(1) = cnf_calc(4) * cnf_calc(5)
           f_sigmoid_der(2) = cnf_in(5) * (1/(cnf_in(pairkind+8)-cnf_in(pairkind+11))) * cnf_calc(5)

        END IF
     END IF
  END FUNCTION f_sigmoid_der

  !----------------------------------------------------------------------------
  ! Window function with and without derivatives
  !----------------------------------------------------------------------------
  FUNCTION f_window(dummy)

     ! win_in(1)   = d0
     ! win_in(2)   = r0
     ! win_in(3)   = c0
     ! win_in(4)   = m
     ! win_in(5)   = n
     !
     ! win_calc(1) = rij
     ! win_calc(2) = rij-d0
     ! win_calc(3) = h0
     ! win_calc(4) = (rij-d0)/r0
     ! win_calc(4) = (-rij+d0)/r0
     ! win_calc(4) = (rij-d0-c0)/(r0-c0)
     ! win_calc(4) = (-rij+d0-c0)/(r0-c0)

     IMPLICIT NONE
     REAL(8) :: f_window, dummy
     REAL(8), DIMENSION(4) :: win_calc

     win_calc(1) = dummy
     win_calc(2) = win_calc(1) - win_in(1) ! win_calc(2) = (rij - d0)

     IF ((win_calc(2).LE.-win_in(3)).OR.(win_calc(2).GE.win_in(3))) THEN
        f_window = 0.d0
     ELSE
        win_calc(3) = 1/(1-(win_in(5)*win_in(2)/ &
                      & ( win_in(4)*(win_in(2)-win_in(3)) )) )
        IF (    (0.d0.LE.win_calc(2)) .AND.(win_calc(2).LE.win_in(2))) THEN
           win_calc(4) = win_calc(2)/win_in(2)
           f_window = (win_calc(3)-1.0d0)* win_calc(4)**win_in(4) + 1.0d0

        ELSE IF ((-win_in(2).LE.win_calc(2)).AND.(win_calc(2).LT. 0.d0)) THEN
           win_calc(4) = -win_calc(2)/win_in(2)
           f_window = (win_calc(3)-1.0d0)* (win_calc(4))**win_in(4) + 1.0d0

        ELSE IF ((win_in(2).LT.win_calc(2)).AND.(win_calc(2).LT. win_in(3))) THEN
           win_calc(4) = (win_calc(2)-win_in(3))/(win_in(2)-win_in(3))
           f_window = win_calc(3)* win_calc(4)**win_in(5)

        ELSE
           win_calc(4) = (-win_calc(2)-win_in(3))/(win_in(2)-win_in(3))
           f_window = win_calc(3)* win_calc(4)**win_in(5)
        END IF
     END IF
  END FUNCTION f_window

  FUNCTION f_window_der(dummy)

     ! win_in(1)   = d0
     ! win_in(2)   = r0
     ! win_in(3)   = c0
     ! win_in(4)   = m
     ! win_in(5)   = n
     !
     ! win_calc(1) = rij
     ! win_calc(2) = rij-d0
     ! win_calc(3) = h0
     ! win_calc(4) = (rij-d0)/r0
     ! win_calc(4) = (-rij+d0)/r0
     ! win_calc(4) = (rij-d0-c0)/(r0-c0)
     ! win_calc(4) = (-rij+d0-c0)/(r0-c0)
     ! win_calc(5) = (h0-1)* ((rij-d0)/r0)**(m-1)
     ! win_calc(5) = (h0-1)* ((-rij+d0)/r0)**(m-1)
     ! win_calc(5) = h0 * (( rij-d0-c0)/(r0-c0))**(n-1)
     ! win_calc(5) = h0 * ((-rij+d0-c0)/(r0-c0))**(n-1)

     ! f_window_der(1) = f_window
     ! f_window_der(2) = derivative of f_window with respect of l_ij

     IMPLICIT NONE
     REAL(8) :: dummy
     REAL(8), DIMENSION(2) :: f_window_der
     REAL(8), DIMENSION(5) :: win_calc

     win_calc(1) = dummy
     win_calc(2) = win_calc(1) - win_in(1) ! win_calc(2) = (rij - d0)

     IF ((win_calc(2).LE.-win_in(3)).OR.(win_calc(2).GE.win_in(3))) THEN
        f_window_der(1:2) = 0.d0
     ELSE
        win_calc(3) = 1/(1-(win_in(5)*win_in(2)/ &
                   & ( win_in(4)*(win_in(2)-win_in(3)) )) )

        IF (    (0.d0.LE.win_calc(2)) .AND.(win_calc(2).LE.win_in(2))) THEN
           win_calc(4) = win_calc(2)/win_in(2)
           win_calc(5) = (win_calc(3)-1.0d0)* win_calc(4)**(win_in(4)-1)

           f_window_der(1) = win_calc(4) * win_calc(5) + 1.0d0
           f_window_der(2) = (1/win_in(2)) * win_in(4) * win_calc(5)

        ELSE IF ((-win_in(2).LE.win_calc(2)).AND.(win_calc(2).LT. 0.d0)) THEN
           win_calc(4) = -win_calc(2)/win_in(2)
           win_calc(5) = (win_calc(3)-1.0d0)* (win_calc(4))**(win_in(4)-1)

           f_window_der(1) = win_calc(4) * win_calc(5) + 1.0d0
           f_window_der(2) = -(1/win_in(2)) * win_in(4) * win_calc(5)

        ELSE IF ((win_in(2).LT.win_calc(2)).AND.(win_calc(2).LT. win_in(3))) THEN
           win_calc(4) = (win_calc(2)-win_in(3))/(win_in(2)-win_in(3))
           win_calc(5) = win_calc(3)* win_calc(4)**(win_in(5)-1)

           f_window_der(1) = win_calc(4) * win_calc(5)
           f_window_der(2) = (1/(win_in(2)-win_in(3))) * win_in(5) * win_calc(5)

        ELSE
           win_calc(4) = (-win_calc(2)-win_in(3))/(win_in(2)-win_in(3))
           win_calc(5) = win_calc(3)* win_calc(4)**(win_in(5)-1)

           f_window_der(1) = win_calc(4) * win_calc(5)
           f_window_der(2) = -(1/(win_in(2)-win_in(3))) * win_in(5) * win_calc(5)

        END IF
     END IF

  END FUNCTION f_window_der

END MODULE META

max_stars_repo_name: kcl-tscm/LoDiS

max_issues_repo_name: kcl-tscm/LoDiS

max_forks_repo_name: kcl-tscm/LoDiS


text:
My names  Jason Kong

Im an avid eaterer of davis restaurants

I tell it how it is


text:
   ! I think this is the same problem as test2007_107.f90
     DOUBLE PRECISION FIELD(5,5)

   ! This appears to be an arror (as reported by gfortran)
     DATA ((FIELD(I,J),I=1,J), J=1,5) / 15 * -1.0D0 /
END

max_stars_repo_name: maurizioabba/rose


text: #redirect Activities and Recreation Center Ballroom