quantum-espresso/PHonon/PH/check_initial_status.f90

657 lines
24 KiB
Fortran

!
! Copyright (C) 2012-2023 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
SUBROUTINE check_initial_status(auxdyn)
!-----------------------------------------------------------------------
!! This routine checks the initial status of the phonon run and sets
!! the variables that control the run, dealing with the image
!! and GRID parallelization features of the phonon code.
!
!! The size of the grid is determined by the following variables:
!! - \(\text{nqs}\): the number of q points;
!! - \(\text{x_q}\): the coordinates of the q points.
!
! - nfs: the number of imaginary frequencies;
! - fiu: which frequencies.
! The flags that control which tensors to calculate.
!
!! In a recover calculation the q grid variables are already known,
!! read from file in \(\texttt{phq_readin}\). In a calculation starting
!! from scratch this routine sets them. The frequencies variables and
!! the tensors flags are read from input.
!! The amount of work to do for each representation of each q
!! point depends on the size of the representation and the
!! order of the small group of q. In a recover calculation
!! these information are on file, when recover=.false. this
!! routine writes the modes and their degeneration on files
!! and calculates the order of the small group of q.
!! The following variables are set:
!! - \(\text{irr_iq}\): for each q point how many irreducible representations;
!! - \(\text{npert_irr_iq}\): how many perturbation per representation and per q;
!! - \(\text{nsymq_iq}\): the order of the small group of q for each q.
!
!! The following variables are set by this routine on the basis of
!! start_irr, last_irr, start_q, last_q, OR of modenum, OR of ifat and
!! atomo:
!! - \(\text{comp_iq}\): TRUE if the q point is calculated in this run;
!! - \(\text{comp_irr_iq}\): TRUE if the representation is calculated in this run;
!! - \(\text{comp_iu}\): TRUE if this frequency is calculated in this run.
! NB: start_iu, last_iu is not yet programmed
!
!! After knowing this info the routine divides the total work among
!! the images (when nimage > 1) INDEPENDENTLY of what has been already
!! calculated and is available on file.
!! Then, when \(\text{recover}\)=TRUE, the routine looks on files for pieces
!! already calculated and sets the array.
!
!! - \(\text{done_irr_iq}\): TRUE if the representation has been already calculated;
!! - \(\text{done_iq}\): TRUE if the q point has been already calculated;
!! - \(\text{done_iu}\): TRUE if already calculated;
!! - \(\text{done_bands_iq}\): TRUE if the bands for the q point are on file.
!
!! If \(\text{recover}\)=FALSE all these array are initialized to FALSE.
!
!! Finally this routine creates a file \(\texttt{fildyn0}\) and writes the
!! q mesh, if this file is not present in the current directory, or if
!! \(\text{recover}\)=FALSE.
!! It also creates a directory for each q inside outdir/_ph#
!! if this directory does not exist and \(\text{lqdir}\)=TRUE.
!
USE io_global, ONLY : stdout
USE control_flags, ONLY : modenum
USE ions_base, ONLY : nat
USE io_files, ONLY : tmp_dir, postfix
USE lsda_mod, ONLY : nspin
USE scf, ONLY : rho
USE disp, ONLY : nqs, x_q, wq, comp_iq, nq1, nq2, nq3, &
done_iq, lgamma_iq
USE qpoint, ONLY : xq
USE control_lr, ONLY : lgamma
USE output, ONLY : fildyn, fildvscf
USE control_ph, ONLY : ldisp, recover, where_rec, rec_code, &
start_q, last_q, current_iq, tmp_dir_ph, &
ext_recover, ext_restart, tmp_dir_phq, lqdir, &
start_irr, last_irr, newgrid, qplot, &
done_zeu, done_start_zstar, done_epsil, &
done_zue, with_ext_images, always_run, trans, &
u_from_file, epsil
USE save_ph, ONLY : tmp_dir_save
USE units_ph, ONLY : iudyn
USE ph_restart, ONLY : check_directory_phsave, check_available_bands,&
allocate_grid_variables, ph_writefile
USE freq_ph, ONLY : current_iu
USE io_rho_xml, ONLY : write_scf
USE mp_images, ONLY : nimage, intra_image_comm
USE io_global, ONLY : ionode, ionode_id
USE io_files, ONLY : prefix , create_directory
USE mp, ONLY : mp_bcast
USE mp_global, ONLY : mp_global_end
USE el_phon, ONLY : elph_mat
USE noncollin_module, ONLY : noncolin, domag
! YAMBO >
USE YAMBO, ONLY : elph_yambo,dvscf_yambo
! YAMBO <
!
USE acfdtest, ONLY : acfdt_is_active, acfdt_num_der
USE dvscf_interpolate, ONLY : ldvscf_interpolate
!
IMPLICIT NONE
!
CHARACTER (LEN=256) :: auxdyn, filename
CHARACTER (LEN=256), EXTERNAL :: trimcheck
CHARACTER (LEN=6 ), EXTERNAL :: int_to_char
LOGICAL :: exst
LOGICAL :: distribute_irrep
INTEGER :: iq, iq_start, ierr
!
tmp_dir=tmp_dir_ph
!
! If this not a recover run, we generate the q mesh. Otherwise at this
! point the code has read the q mesh from the files contained in
! prefix.phsave
!
IF (.NOT.recover) THEN
!
! recover file not found or not looked for
!
current_iu=1
current_iq=1
IF (ldisp) THEN
!
! ... Calculate the q-points for the dispersion
!
IF(elph_mat) then
CALL q_points_wannier()
ELSE
IF (.NOT. qplot) CALL q_points()
END IF
!
! YAMBO >
ELSE IF (.NOT.elph_yambo .AND. .NOT. dvscf_yambo) then
! YAMBO <
!
nqs = 1
last_q = 1
ALLOCATE(x_q(3,1))
ALLOCATE(wq(1))
ALLOCATE(lgamma_iq(1))
x_q(:,1)=xq(:)
wq(1)=1.0d0
lgamma_iq(1)=lgamma
!
END IF
!
! Save the mesh of q and the control flags on file
!
CALL ph_writefile('init',0,0,ierr)
!
! Initialize the representations and write them on file.
!
IF (trans .OR. epsil .OR. ldvscf_interpolate.or.elph_mat) THEN
CALL init_representations()
ELSE
u_from_file = .TRUE.
ENDIF
!
IF ((start_irr==0).AND.(last_irr==0)) THEN
where_rec='init_rep..'
rec_code=-50
current_iq=1
current_iu=1
CALL ph_writefile('status_ph',current_iq,0,ierr)
CALL clean_pw(.FALSE.)
CALL close_files(.FALSE.)
CALL mp_global_end()
STOP
ENDIF
ENDIF
IF (last_q<1.or.last_q>nqs) last_q=nqs
IF (start_q<1.or.start_q>last_q) call errore('check_initial_status', &
'wrong start_q',1)
!
! now we allocate the variables needed to describe the grid
!
CALL allocate_grid_variables()
!
! This routine assumes that the modes are on file, either written by
! init_representation or written by a previous run. It takes care
! of dealing with start_irr, last_irr flags and ifat or modenum
! restricted computation, moreover it sets the size of each
! representation and the size of the small group of q for each point.
!
CALL initialize_grid_variables()
!
! If there are more than one image, divide the work among the images
!
! If writing dvscf file, distribute only q points to images.
! Otherwise, distribute both q points and irreps to images.
!
distribute_irrep = .TRUE.
IF (fildvscf /= ' ') THEN
WRITE(stdout, '(5x,a)')
WRITE(stdout, '(5x,a)') 'Saving dvscf to file. Distribute only q points, &
&not irreducible representations.'
distribute_irrep = .FALSE.
ENDIF
!
IF (nimage > 1 .AND. .NOT. with_ext_images) CALL image_q_irr(distribute_irrep)
!
IF (recover) THEN
!
! ... Checking the status of the calculation
!
! sets which q point and representations have been already calculated
!
CALL check_directory_phsave()
!
! If a recover or a restart file exists the first q point is the current one.
!
IF ((.NOT.lgamma_iq(current_iq).OR. newgrid).AND.lqdir) THEN
tmp_dir_phq= trimcheck ( TRIM (tmp_dir_ph) // TRIM(prefix) // &
& '.q_' // int_to_char(current_iq) )
tmp_dir=tmp_dir_phq
CALL check_restart_recover(ext_recover, ext_restart)
tmp_dir=tmp_dir_ph
ELSE
CALL check_restart_recover(ext_recover, ext_restart)
ENDIF
IF (.NOT.ext_recover.AND..NOT.ext_restart) THEN
current_iq=start_q
ELSE
!
! Check that the representations from start_q to current_iq have been done
!
DO iq=start_q, current_iq-1
IF (comp_iq(iq) .AND. .NOT.done_iq(iq)) &
CALL errore('check_initial_status',&
& 'recover file found, change in start_q not allowed',1)
comp_iq(iq)=.FALSE.
ENDDO
ENDIF
iq_start=current_iq
!
! check which bands are available and set the array done_bands
!
CALL check_available_bands()
!
! write the information on output
IF (iq_start<=last_q.AND.iq_start>0) THEN
WRITE(stdout, &
'(5x,i4," /",i4," q-points for this run, from", i3,&
& " to", i3,":")') last_q-iq_start+1, nqs, iq_start, last_q
WRITE(stdout, '(5x," N xq(1) xq(2) xq(3) " )')
DO iq = 1, nqs
WRITE(stdout, '(5x,i3, 3f14.9,l6)') iq, x_q(1,iq), x_q(2,iq), &
x_q(3,iq)
END DO
WRITE(stdout, *)
ELSEIF (iq_start>last_q) THEN
WRITE(stdout, &
'(5x,"Starting q",i4," larger than total number of q points", i4, &
& " or of last q ", i3)') iq_start, nqs, last_q
ELSEIF (iq_start<0) THEN
CALL errore('check_initial_status','wrong iq_start',1)
ENDIF
ELSE
done_zeu=.FALSE.
done_start_zstar=.FALSE.
done_epsil=.FALSE.
done_zue=.FALSE.
ENDIF
!
! Create a new directory where the ph variables are saved and copy
! the charge density there.
!********************** ACFDT TEST ***************
IF (acfdt_is_active) THEN
! ACFDT -test always write rho on file
IF (acfdt_num_der) THEN
CALL write_scf( rho, nspin )
ELSE
IF ((ldisp.OR..NOT.lgamma.OR.modenum/=0).AND.(.NOT.lqdir)) &
CALL write_scf( rho, nspin )
ENDIF
ELSE
! this is the standard treatment
IF ( ( ( ldisp.OR..NOT.lgamma .OR. modenum/=0 ) .AND. (.NOT.lqdir) ) &
.OR. newgrid .OR. always_run .OR. (noncolin.AND.domag) ) CALL write_scf( rho, nspin )
ENDIF
!*********************** END OF ACFDT TEST *****************
!
! Write the file fildyn0 with the mesh of q points. This file is used
! by postprocessing programs such as q2r.x and we write it again if
! it is not found in the running directory.
!
filename=TRIM(fildyn)//'0'
ierr=0
IF (ionode.and..NOT.elph_mat) THEN
INQUIRE (FILE = TRIM(filename), EXIST = exst)
ierr=0
IF ((.NOT. exst .OR. .NOT. recover).AND.ldisp) THEN
iudyn=26
OPEN (unit=iudyn, file=TRIM(filename), status='unknown', iostat=ierr)
IF ( ierr == 0 ) THEN
WRITE (iudyn, '(3i4)' ) nq1, nq2, nq3
WRITE (iudyn, '( i4)' ) nqs
DO iq = 1, nqs
WRITE (iudyn, '(3e24.15)') x_q(1,iq), x_q(2,iq), x_q(3,iq)
END DO
CLOSE (unit=iudyn)
ENDIF
ENDIF
END IF
CALL mp_bcast(ierr, ionode_id, intra_image_comm)
CALL errore ('check_initial_status','cannot open file ' // TRIM(filename),&
abs(ierr))
!
! The following commands deal with the flag lqdir=.true. In this case
! each q point works on a different directory. We create the directories
! if they do not exist and copy the self consistent charge density
! there.
!
DO iq = 1,nqs
IF (.NOT.comp_iq(iq)) CYCLE
lgamma = lgamma_iq(iq)
!
! ... each q /= gamma works on a different directory. We create them
! here and copy the charge density inside
!
IF ((.NOT.lgamma.OR. newgrid .OR. (qplot .AND. iq /=1)) .AND.lqdir) THEN
tmp_dir_phq= trimcheck ( TRIM (tmp_dir_ph) // TRIM(prefix) // &
& '.q_' // int_to_char(iq) )
#if defined(__HDF5)
filename=TRIM(tmp_dir_phq)//TRIM(prefix)//postfix//'charge-density.hdf5'
#else
filename=TRIM(tmp_dir_phq)//TRIM(prefix)//postfix//'charge-density.dat'
#endif
IF (ionode) inquire (file =TRIM(filename), exist = exst)
!
CALL mp_bcast( exst, ionode_id, intra_image_comm )
!
IF (.NOT. exst) THEN
CALL create_directory( tmp_dir_phq )
tmp_dir=tmp_dir_phq
CALL write_scf( rho, nspin )
tmp_dir=tmp_dir_save
ENDIF
ENDIF
ENDDO
!
auxdyn = fildyn
RETURN
END SUBROUTINE check_initial_status
SUBROUTINE image_q_irr(distribute_irrep)
!
!! This routine is an example of the load balancing among images.
!! It decides which image makes which q and which irreducible representation.
!! The algorithm at the moment is straightforward. Possibly better
!! methods could be found.
!! It receives as input:
!! - \(\text{nsym}\): the dimension of the point group;
!! - \(\text{nsymq_iq}\): the dimension of the small group of q for each q;
!! - \(\text{irr_iq}\): the number of irreps for each q;
!! - \(\text{npert_irr_iq}\): for each q and each irrep its dimension.
!! It provides as output the two arrays:
!! - \(\text{comp_iq}\): if this q has to be calculated by the present image;
!! - \(\text{comp_irr_iq}\): for each q the array to be copied into comp_irr.
USE ions_base, ONLY : nat
USE disp, ONLY : comp_iq, nqs, nq1, nq2, nq3
USE grid_irr_iq, ONLY : irr_iq, npert_irr_iq, comp_irr_iq, nsymq_iq
USE control_ph, ONLY : start_q, last_q
USE io_global, ONLY : stdout
USE mp_images, ONLY : nimage, my_image_id
USE symm_base, ONLY : nsym
USE modes, ONLY : nmodes
IMPLICIT NONE
LOGICAL, INTENT(IN) :: distribute_irrep
!! If TRUE, distribute both irreps and q-points (default behavior).
!! If FALSE, distribute only q-points. Used when writing dvscf file.
!
INTEGER :: total_work, & ! total amount of work to do
total_nrapp, & ! total number of representations
total_nq, & ! total number of q points
work_per_image ! approximate minimum work per image
INTEGER, ALLOCATABLE :: image_iq(:,:), work(:)
INTEGER :: iq, irr, image, work_so_far, actual_diff, diff_for_next
CHARACTER(LEN=256) :: string
CHARACTER(LEN=6), EXTERNAL :: int_to_char
ALLOCATE (image_iq(0:3*nat,nqs))
ALLOCATE (work(0:nimage-1))
IF (distribute_irrep) THEN
total_work=0
total_nrapp=0
DO iq = start_q, last_q
DO irr = 1, irr_iq(iq)
IF (comp_irr_iq(irr,iq)) THEN
total_work = total_work + npert_irr_iq(irr, iq) * nsym / nsymq_iq(iq)
IF (irr==1) total_work = total_work + nsym / nsymq_iq(iq)
total_nrapp = total_nrapp + 1
ENDIF
END DO
END DO
IF (nimage > total_nrapp) &
CALL errore('image_q_irr','some images have no rapp', 1)
work_per_image = total_work / nimage
!
! If nimage=total_nrapp we put one representation per image
! No load balancing is possible. Otherwise we try to minimize the number of
! different q per image doing all representations of a given q until
! the work becomes too large.
! The initialization is done by the image with the first representation of
! each q point.
!
image=0
work=0
work_so_far=0
DO iq = start_q, last_q
DO irr = 1, irr_iq(iq)
IF (comp_irr_iq(irr,iq)) THEN
image_iq(irr,iq) = image
work(image)=work(image) + npert_irr_iq(irr, iq) * nsym / nsymq_iq(iq)
work_so_far=work_so_far + npert_irr_iq(irr, iq) * nsym / nsymq_iq(iq)
IF (irr==1) THEN
image_iq(0,iq)=image
work(image)=work(image) + nsym / nsymq_iq(iq)
work_so_far=work_so_far + nsym / nsymq_iq(iq)
ENDIF
!
! The logic is the following. We know how much work the current image
! has already accumulated and we calculate how far it is from the target.
! Note that actual_diff is a positive number in the usual case in which
! we are below the target. Then we calculate the work that the current
! image would do if we would give it the next representation. If the work is
! still below the target, diff_for_next is negative and we give the
! representation to the current image. If the work is above the target,
! we give it to the current image only if its distance from the target
! is less than actual_diff.
!
actual_diff=-work(image)+work_per_image
IF (irr<irr_iq(iq)) THEN
diff_for_next= work(image)+npert_irr_iq(irr+1, iq)*nsym/nsymq_iq(iq) &
- work_per_image
ELSEIF (irr==irr_iq(iq).and.iq<last_q) THEN
diff_for_next= work(image)+npert_irr_iq(1, iq+1)* &
nsym/nsymq_iq(iq+1) + nsym/nsymq_iq(iq+1)-work_per_image
ELSE
diff_for_next=0
ENDIF
IF ((nimage==total_nrapp.OR.diff_for_next>actual_diff).AND. &
(image < nimage-1)) THEN
work_per_image= (total_work-work_so_far) / (nimage-image-1)
image=image+1
ENDIF
ENDIF
ENDDO
ENDDO
!
ELSE ! .NOT. distribute_irrep
!
! We distribute only q-points.
! The work for each q-point is (nmodes + 1) * nsym / nsymq_iq(iq).
!
total_work = 0
total_nq = 0
DO iq = start_q, last_q
total_work = total_work + (nmodes + 1) * nsym / nsymq_iq(iq)
total_nq = total_nq + 1
ENDDO
IF (nimage > total_nq) &
CALL errore('image_q_irr','some images have no rapp', 1)
work_per_image = total_work / nimage
!
! If nimage=total_nq we put one representation per image
! No load balancing is possible. Otherwise we try to minimize the number of
! different q per image doing all representations of a given q until
! the work becomes too large.
! The initialization is done by the image with the first representation of
! each q point.
!
image = 0
work = 0
work_so_far =0
DO iq = start_q, last_q
!
! Assign iq to this image
work(image) = work(image) + (nmodes + 1) * nsym / nsymq_iq(iq)
work_so_far = work_so_far + (nmodes + 1) * nsym / nsymq_iq(iq)
!
! Assign all irrs of iq to this image
DO irr = 1, irr_iq(iq)
IF (comp_irr_iq(irr, iq)) THEN
image_iq(irr, iq) = image
ENDIF
ENDDO
image_iq(0, iq) = image
!
! See the comment above for an explanation of the logic.
!
actual_diff = work_per_image - work(image)
IF (iq < last_q) THEN
diff_for_next = work(image) + (nmodes + 1) * nsym / nsymq_iq(iq + 1) &
- work_per_image
ELSE
diff_for_next = 0
ENDIF
IF ((nimage == total_nq .OR. diff_for_next > actual_diff) &
.AND. (image < nimage - 1)) THEN
work_per_image = (total_work - work_so_far) / (nimage - image - 1)
image = image + 1
ENDIF
ENDDO ! iq
!
ENDIF ! distribute_irrep
!
! Here we actually distribute the work. This image makes only
! the representations calculated before.
!
DO iq = start_q, last_q
DO irr = 0, irr_iq(iq)
IF (image_iq(irr,iq)/=my_image_id ) THEN
comp_irr_iq(irr,iq)=.FALSE.
ENDIF
ENDDO
ENDDO
comp_iq = .FALSE.
DO iq = start_q, last_q
DO irr = 0, irr_iq(iq)
IF (comp_irr_iq(irr,iq).AND..NOT.comp_iq(iq)) THEN
comp_iq(iq)=.TRUE.
ENDIF
ENDDO
ENDDO
WRITE(stdout, &
'(/,5x," Image parallelization. There are", i3,&
& " images", " and ", i5, " representations")') nimage, &
total_nrapp
WRITE(stdout, &
'(5x," The estimated total work is ", i5,&
& " self-consistent (scf) runs")') total_work
WRITE(stdout, '(5x," I am image number ",i5," and my work is about",i5, &
& " scf runs. I calculate: ")') &
my_image_id, work(my_image_id)
DO iq = 1, nqs
IF (comp_iq(iq)) THEN
WRITE(stdout, '(5x," q point number ", i5, ", representations:")') iq
string=' '
DO irr=0, irr_iq(iq)
IF (comp_irr_iq(irr, iq)) &
string=TRIM(string) // " " // TRIM(int_to_char(irr))
ENDDO
WRITE(stdout,'(6x,A)') TRIM(string)
ENDIF
ENDDO
DEALLOCATE(image_iq)
DEALLOCATE(work)
RETURN
END SUBROUTINE image_q_irr
SUBROUTINE collect_grid_files()
!
!! This subroutine collects all the xml files contained in different
!! directories and created by the diffent images in the phsave directory
!! of the image 0.
!
USE io_files, ONLY : tmp_dir, prefix
USE control_ph, ONLY : tmp_dir_ph
USE save_ph, ONLY : tmp_dir_save
USE disp, ONLY : nqs, lgamma_iq
USE grid_irr_iq, ONLY : comp_irr_iq, irr_iq
USE control_ph, ONLY : ldisp, epsil, zue, zeu
USE klist, ONLY : lgauss, ltetra
USE el_phon, ONLY : elph
USE clib_wrappers, ONLY : f_copy
USE mp, ONLY : mp_barrier
USE mp_images, ONLY : my_image_id, nimage, intra_image_comm
USE io_global, ONLY : stdout, ionode
IMPLICIT NONE
INTEGER :: iq, irr, ios
LOGICAL :: exst
CHARACTER(LEN=256) :: file_input, file_output
CHARACTER(LEN=6), EXTERNAL :: int_to_char
CHARACTER(LEN=8) :: phpostfix
#if defined (_WIN32)
#if defined (__PGI)
phpostfix='.phsave\\'
#else
phpostfix='.phsave\'
#endif
#else
phpostfix='.phsave/'
#endif
CALL mp_barrier(intra_image_comm)
IF (nimage == 1) RETURN
IF (my_image_id==0) RETURN
DO iq=1,nqs
DO irr=0, irr_iq(iq)
IF (comp_irr_iq(irr,iq).and.ionode) THEN
file_input=TRIM( tmp_dir_ph ) // &
& TRIM( prefix ) // phpostfix // 'dynmat.' &
& // TRIM(int_to_char(iq))&
& // '.' // TRIM(int_to_char(irr)) // '.xml'
file_output=TRIM( tmp_dir_save ) // '/_ph0/' &
& // TRIM( prefix ) // phpostfix // 'dynmat.' &
& // TRIM(int_to_char(iq)) &
& // '.' // TRIM(int_to_char(irr)) // '.xml'
INQUIRE (FILE = TRIM(file_input), EXIST = exst)
IF (exst) ios = f_copy(file_input, file_output)
IF ( elph .AND. irr>0 ) THEN
file_input=TRIM( tmp_dir_ph ) // &
& TRIM( prefix ) // phpostfix // 'elph.' &
& // TRIM(int_to_char(iq))&
& // '.' // TRIM(int_to_char(irr)) // '.xml'
file_output=TRIM( tmp_dir_save ) // '/_ph0/' // &
& TRIM( prefix ) // phpostfix // 'elph.' &
& // TRIM(int_to_char(iq)) &
& // '.' // TRIM(int_to_char(irr)) // '.xml'
INQUIRE (FILE = TRIM(file_input), EXIST = exst)
IF (exst) ios = f_copy(file_input, file_output)
ENDIF
ENDIF
ENDDO
IF ((ldisp.AND..NOT. (lgauss .OR. ltetra)).OR.(epsil.OR.zeu.OR.zue)) THEN
IF (lgamma_iq(iq).AND.comp_irr_iq(0,iq).AND.ionode) THEN
file_input=TRIM( tmp_dir_ph ) // &
TRIM( prefix ) // phpostfix // 'tensors.xml'
file_output=TRIM( tmp_dir_save ) // '/_ph0/' &
// TRIM( prefix ) // phpostfix // 'tensors.xml'
INQUIRE (FILE = TRIM(file_input), EXIST = exst)
IF (exst) ios = f_copy(file_input, file_output)
ENDIF
ENDIF
ENDDO
RETURN
END SUBROUTINE collect_grid_files