SUBROUTINE aph3d1 (rmat, eye11, w, vecnew, tstore, vecold, xsq, xk, nidm) ! ! P U R P O S E O F S U B R O U T I N E ! This is the Main program for the APH3D1 propagator. ! I N P U T A R G U M E N T S ! O U T P U T A R G U M E N T S ! <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> USE FileUnits_Module USE nstate_Module USE Narran_Module USE parms_Module USE Masses_Module USE TotalEng_Module USE Regins_Module USE Title1_Module USE region_MODULE USE Storage_Module USE convrsns_MODULE USE aphchl_MODULE USE VFunc_Module USE NumNuj_Module USE chltot_MODULE USE engpro_MODULE USE APHchl_Module, ONLY: Neigmin_Sum IMPLICIT NONE REAL(Kind=WP_Kind) version, eredev REAL(Kind=WP_Kind) rmat(nstate,nstate), eye11(nstate,nstate), w(nstate,nstate) REAL(Kind=WP_Kind) vecnew(nstate,nstate), tstore(nstate,nstate), vecold(nstate,nstate) REAL(Kind=WP_Kind) xsq(nstate), xk(nstate) INTEGER ithcall, ithsub, nidm, njacobi,i CHARACTER(LEN=10) hour CHARACTER(LEN=8) today CHARACTER(LEN=5) curzone INTEGER dtvalues(8) LOGICAL iaph, idelve, ijacobi, little, medium, full EXTERNAL popt, date, time, JacBasis, vsets, aphprop, delprop, jacprop DATA iaph/.false./, idelve/.false./, ijacobi/.false./ !------------------------------------------------------------------ ! Determine printing options. !------------------------------------------------------------------ DATA little/.false./,medium/.false./,full/.false./, ithcall/0/,ithsub/0/,version/3.3/ CALL popt('aph3d1 ',little,medium,full,ithcall,ithsub) !------------------------------------------------------------------ ! Call system routines to determine the date and time. !------------------------------------------------------------------ CALL dateh(today) CALL timeh(hour) !------------------------------------------------------------------ ! WRITE the titles for this particular run. !------------------------------------------------------------------ WRITE(Out_Unit,1131) WRITE(Out_Unit,*)'These results were made with ', 'APH3D1 version ',version !------------------------------------------------------------- ! WRITE the time and date to unit 6. !------------------------------------------------------------- CALL Print_Date_And_Time(Today, Hour, CurZone, DtValues, Out_Unit) WRITE(Out_Unit,1132)Title1 WRITE(Out_Unit,1132)Title2 WRITE(Out_Unit,1132)Title3 !------------------------------------------------------------------ ! The maximum energy will determine the wavelength ! of the most oscillatory state which will be used to set the ! step size in the log-derivative propagation. !------------------------------------------------------------------ engmax=emin+(numengs-1)*ejump !------------------------------------------------------------------ ! If irdindep=true the energy independent information has already ! been stored on the disk and hence iwrindep=false. ! If irdindep=false we usually WRITE the energy independent stuff ! to disk, but if only doing one energy both irdindep ! and iwrindep can be false. !------------------------------------------------------------------ IF(irdindep) iwrindep=.false. ! ***************************** ! * loop for each energy. * ! ***************************** kenergy = 1 !------------------------------------------------------------------ ! When irdindep=true the energy independent information will be ! read from disk. This implies that a previous run has been made with ! iwrindep=true, irdindep=false. !------------------------------------------------------------------ IF(kenergy>1)THEN irdindep=.true. iwrindep=.false. OPEN(Unit=Basis_Bin_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'BinOut/Basis.bin',form='unformatted',status='old') REWIND frst_unit ELSE IF(iwrindep)THEN OPEN(Unit=Basis_Bin_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'BinOut/Basis.bin',form='unformatted', status='unknown') WRITE(frst_unit)usys2 ELSEIF(irdindep)THEN OPEN(Unit=Basis_Bin_Unit,File=OutDIR(1:LEN(TRIM(OutDIR)))//'BinOut/Basis.bin',form='unformatted',status='old') ENDIF ENDIF !------------------------------------------------------------------ ! Determine the Total collision energy. !------------------------------------------------------------------ Etot=emin+(kenergy-1)*ejump eredev=Etot*autoev WRITE(Out_Unit,184)Etot,eredev 184 FORMAT(//,'energy (Ha) =',f10.5,5x,'energy (ev) =',f10.5) !-------------------------------------------------------------------- ! Obtain the asymptotic Jacobi basis functions. !-------------------------------------------------------------------- WRITE(Out_Unit,*)' G E N E R A T E J A C O B I B A S I S' ! CALL JacBasis (nchanl, vecold, tstore, xsq, xk, finish, chinuj, nnuj, njacobi) WRITE(Out_Unit,*)'finished jacobi basis' !-------------------------------------------------------------------- ! Set the initial R-matrix to essentially a zero matrix. ! The diagonal elements are set to a very small number to ! avoid overflow problems in inverting the R-matrix to ! obtain the Log-Derivative matrix. !-------------------------------------------------------------------- naph=Neigmin_Sum rmat=0.0d0 DO i=1,naph rmat(i,i)=1d-30 ENDDO IF(naph>nstate)THEN Write(Msg_Unit,*)'Error. naph greater than nstate' STOP 'aph3d1' ENDIF !-------------------------------------------------------------------- ! Propagate the coupled equations in the APH region from start ! to endaph if this region has nonzero lenght. !-------------------------------------------------------------------- WRITE(Out_Unit,*)start,endaph IF(start<=endaph)THEN iregion='aph ' iaph=.true. WRITE(Out_Unit,*)' A P H P R O P A G A T I O N' CALL aphprop (naph,rmat,eye11,w,vecnew, tstore,vecold,xsq,start,endaph,njacobi,xk,nidm) WRITE(Out_Unit,*)'END APH PROPAGATION' ENDIF !-------------------------------------------------------------------- ! Propagate the coupled equations in the Delves region from ! endaph to enddelve if this region has nonzero length. !-------------------------------------------------------------------- WRITE(Out_Unit,*)endaph,enddelve IF(endaph