module energylevs ! this is the module which stores the energy level quantum numbers ! and energies implicit none save DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: enlevs integer, allocatable, dimension(:) :: ij, isym4, inn integer, allocatable, dimension(:,:) :: iV, KJ integer nlevs end module energylevs !------------------------------ program spect4 ! this program is designed to handle the output from programs dipole ! (of the TRIATOM suite) and dipole3 (of the DVR3D suite) ! for several runs and turn it into data sorted on the transition ! frequency wif. From this, and input from triatom/dvr3d and rotlev runs, ! the partition functions and integrated absorption coefficients at ! the required temperatures. ! the program structure consists essentially of four main parts: ! 1) the main program, spect, determines the amount of memory needed ! for the job in hand, and calls subroutine gtmain, a machine ! dependent routine, which allocates the memory needed. spect also ! calls the other three main working routines. ! 2) subroutine sortsp. this routine sorts out data from program ! dipole on ascending transition frequency and may print this out ! if required. ! the main input file is "itra", and contains the output from the ! necessary dipole runs. by default itra is attached to stream 13. ! itra is read in the main program to determine the size of the ! problem. the main program also creates a scratch file, item (16), ! to cut down the amount of data sorted in sortsp. ! subroutine sortsp creates a scratch file ispe for calculations of ! spectral intensities, if required. if this file has been saved from ! a previous run, sortsp may be skipped. ispe defaults to stream 15. ! 3) subroutine pfcalc. this routine calculates the partition function ! from energy levels calculated by programs triatom and rotlevd at the ! temperature required for each run. ! data necessary to calculate the partition function is given on ! file "ilev". this is defaulted to stream 14. ! the first eigenenergy on ilev must be the ground state for j=0, ! (q=0). this is the zero energy of the problem. ilev is read in the ! main program to determine the size of the problem and in subroutine ! pfcal! to access the data. ! 4) subroutine spectm. this routine calculates boltzman weighted ! intensities at the required temperatures and in the frequency ranges ! determined by the input parameters given by input data lines 4-n. ! spectm reads the sorted data from sortsp off scratch file ispe. ! if required it writes out data for a suitable graphics package to ! data stream iplot (default 20). one option is to produce a spectrum ! of gaussian profiles by calling subroutine profil. users may wish to ! substitute this routine with other profile types. ! lines can also be written to a linelist file, stream ilist (default ! 36), if zlist is .true. ! input data is supplied to the program on stream 5, as follows: ! 1) namelist prt: four logical parameters: ! zout(false) this should be true if the sorted ! line strengths are to be written ! to the lineprinter. zout is set to ! true automatically if zspe is false. ! zsort(true) if false subroutine sortsp is skipped. ! zspe (true) if false the program stops after sortsp. ! units of ispe are atomic units. ! zpfun(true) calculates the partition function ! from energy levels supplied from ! DVR3DRJZ and ROTLEV3/3B. ! if zpfun false, the partition function ! is set to q read in below. ! ZBAND(false) Set this to true if a specific band is required ! Must then provide upper and lower ! vibrational quantum numbers. ! ! the default values of itra(13), ilev(14), ispe(15), ! and item(16) may also be reset on this namelist ! if required. ! gz(0.0) absolute energy of the ground state ! in cm-1 (not used if zpfun=.true.). ! If zsort=.true. the range of data to be sorted may be set by setting ! any or all of wsmin, wsmax, emin, emax, jmax, smin. ! emin(-1.d27) the minimum value of e" for which data ! is to be written out. ! emax(+1.d27) the maximum value of e" for which data ! is to be written out. if emax is not ! specified all values of e" above emin ! are considered. ! jmax(-1) the maximum value of j" for which ! transition frequencies, et! are to be ! written out in subroutine sortspe. ! if jmax is not specified, all j" values ! are considered. ! smin(0.0d0) The minimum value of the square of the ! transition dipole (D**2). ! ! smax The maximum value of the square of the ! transition dipole (D**2) ! ! wsmax The maximum frequency of a transition (cm-1). ! ! wsmin The minimum frequency of a transition (cm-1) ! judicious use of the parameters on this namelist can enable certain types ! of transitions to be separated - e.g. fundamental bands, hot bands etc. ! 2) a title, format 9a8. ! 3) ge, go. format 2f10.0 nuclear degeneracy factors for homonuclear ! diatomic molecules for the j even and j odd states. ! defaults of 1.0d0 supplied if left totally blank. ! if the molecule does not contain a homonuclear diatomic ! part, ge and go are ignored. ! 4) temp, xmin, wmin, wmax, dwl, q. format 6f10.0 ! temp temperature in k; ! xmin the lowest relative intensity to be printed out; ! wmin and wmax define the transition frequency range in cm-1. ! dwl the width of the gaussian profile for the lines (see below). ! q is the partition function (only used if ZPFUN=F). ! use of default (q=1.0) means the absolute absortpion ! coefficients are not correct. ! 5) namelist spe. this namelist is read in subroutine spectm. ! emin1, emax1, emin2, emax2, tinte and jmax control the print out ! of spectral intensities. ! tinte set the value below which the intensities are printed out. ! zemit (.false.) gives emissivities if .true. ! in addition, there is a parameter zplot(false) which creates an ! ascii file suitable for use with a graphics package. this can ! consist of stick spectra or a gaussian broadened spectrum if ! zprof (.false.) is .true. at npoint points, it is in: ! 1) frequency in wavenumbers against intensity if zfreq (.true.) ! is .true. ! 2) wavelength in microns against intensity if zfreq .false. ! 3) either of the above against einstein a-coefficient times ! spin weighting in zeinst (.false.) is .true. ! zeinst (.true.) sets zemit to .true. ! in addition if zene is .true the energy levels and the rotational ! levels are printed in the output file iplot. ! if zprof is .true., stream idat takes data from spectm to profil. ! The program makes use of the dpsort freeware routines, found at ! http://netlib2.cs.utk.edu/ ! ! Updated to f90 to use dynamic memory allocation and to preselect transitions ! by GJH & JT 2001. use energylevs implicit double precision(a-h,o-y), logical(z) ! DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: enlevs ! integer, allocatable, dimension(:) :: ij, isym4, inn ! integer, allocatable, dimension(:,:) :: VV, KJ character(len=8) title(9) namelist/prt/ zout, zsort, zspe, zpfun, itra, ilev, ispe, item, & wsmax,wsmin, emin, emax, jmax, smin, gz, zband common/logic/ zout, zsort, zspe, zpfun, zembed common/base/ ibase1,ibase2 common/timing/itime0 ! common/energylevs/ enlevs,ij,isym4,inn, vv,kj, nenlevs ! DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: ee1, ee2, ss data autocm/ 2.19474624d+05/, autode/ 2.5417662d0/ ! detosc converts from s(f-i) in debye**2 to seconds**(-1) data detosc/ 3.136186d-07/ nlevs=163491 !(hdo) this is the total number of energy levels *******************************our ! set time zero for timer call SYSTEM_CLOCK(itime0,irate2,imax2) ! preset namelist defaults zout= .false. zsort= .true. zspe= .true. zpfun= .true. zband= .false. itra= 13 ilev= 14 ispe= 15 item= 16 wsmin= 0.0d0 wsmax= 1.0d6 smin=0.0d0 jmax=500 emin= -1.0d27 emax= +1.0d27 gz=0.0d0 ! NOTE IDIA IS HARD CODED HERE idia=-2 write(6,*) "Program SPECTRA (version of January 2004)" ! read in logical variables and title open(unit=5,file='spec.job',form='formatted') read(5,prt) read(5,100) title 100 format(9a8) write(6,202) title 202 format(/5x,9a8/) if (zband) read(5,*) ivu1, ivu2, ivu3, ivl1, ivl2, ivl3 ! allocate memory of the energy levels allocate(enlevs(nlevs),ij(nlevs),isym4(nlevs),inn(nlevs)) allocate(kj(nlevs,3), iv(nlevs,3)) call read_energies ! determine energy zero from levels file if available if (zpfun) then gz=enlevs(1) endif write(6,203) gz 203 format(5x,'Ground state energy set to',f14.7,' cm-1') ! determine number of transitions to be processed and create item maxr=0 nr= 0 nmax1= 0 nmax2= 0 open(unit=ispe,form='unformatted') if(zsort) then write(6,1000) wsmin,wsmax,emin,emax,smin,jmax 1000 format(/5x,'Sorting requested, list preselected using:', & /5x,'Minimum frequency, WSMIN =',d12.5,' cm-1', & /5x,'Maximum frequency, WSMAX =',d12.5,' cm-1', & /5x,'Minimum energy, EMIN =',d12.5,' cm-1', & /5x,'Maximum energy, EMAX =',d12.5,' cm-1', & /5x,'Minimum linestrength, SMIN =',d12.5,' D**2', & /5x,'Maximum rotations JMAX =',i5) wsmi=wsmin/autocm wsma=wsmax/autocm emi=(emin+GZ)/autocm ema=(emax+GZ)/autocm smi=smin/autode**2 ! begining of 1st select phase ! ! This phase reads in the data line by line and calculates ! S, frequncy and other parameters then if criterior are met ! outputs to item. ! open(unit=itra,form='formatted') open(unit=item,form='unformatted') rewind itra maxr=0 nr=0 10 read(itra,102,end=90) index1,index2,AA maxr=maxr+1 EE1=enlevs(index1) EE2=enlevs(index2) if(zband) then if(iv(index1,1) .eq. ivu1 .and. iv(index1,2) .eq. ivu2 .and. & iv(index1,3) .eq. ivu3 .and. iv(index2,1) .eq. ivl1 .and. & iv(index2,2) .eq. ivl2 .and. iv(index2,3) .eq. ivl3) goto 15 if(iv(index1,1) .eq. ivl1 .and. iv(index1,2) .eq. ivl2 .and. & iv(index1,3) .eq. ivl3 .and. iv(index2,1) .eq. ivu1 .and. & iv(index2,2) .eq. ivu2 .and. iv(index2,3) .eq. ivu3) goto 15 goto 10 endif 15 if (ee1.ge.emin .and. ee1.le.emax) then if (ee2.ge.emin .and. ee2.le.emax) then w= ee2 - ee1 if(w.eq.0.000000) w=0.0000001 !!!!!!!!!!!!!!! fix ZERO Boris July 2009 if (abs(w).ge.wsmin .and. abs(w).le.wsmax) then !! write(6,*) 'isym4(index1)=',isym4(index1), 'ij(index1)=', ij(index1) !! Boris++++++++++++++++++++++++++++++++++++++++++ !! write(6,*) 'isym4(index2)=',isym4(index2), 'ij(index2)=', ij(index2) !! ++++++++++++++++++++++++++++++++++++++++++ call isym_reverse(isym4(index1), ij(index1), ipar1, kmin1) call isym_reverse(isym4(index2), ij(index2), ipar2, kmin2) if (w.ge.0.0d0) then ss=(dble(2*ij(index2) + 1)*aa)/abs(w*w*w*detosc) else ss=-(dble(2*ij(index1) + 1)*aa)/(w*w*w*detosc) endif !-- if (ss.ge.smin) then ss=ss/autode**2 ee1=ee1/autocm ee2=ee2/autocm w=w/autocm nr=nr+1 if (w.ge.0.0d0) then write(item) ipar1,ij(index2),kmin2,inn(index2), & ij(index1),kmin1,inn(index1),index2,index1, & ee2,ee1,w,ss,1,1 else write(item) ipar1,ij(index1),kmin1,inn(index1),& ij(index2),kmin2,inn(index2),index1,index2, & ee1,ee2,-w,ss,1,1 endif endif endif endif endif goto 10 90 close(itra) write(6,1010) maxr,itra,nr,item 1010 format(/,i10,' transition records read from unit ITRA =',i3, & /,i10,' selected and written to unit ITEM =',I3) if (nr.le.0) then write(6,900) 900 format(/'No lines selected: stop') stop endif else rewind ispe read(ispe) 13 read(ispe,end=93) nr= nr+1 goto 13 93 continue write(6,1020) nr,ispe 1020 format(/5x,'List already sorted:', & /i10,' transition records found on unit ISPE =',i3) endif write (6,*) "call spmain" call spmain(ilev,ispe,nr,nmax1,nmax2,item,idia,gz) deallocate(enlevs,ij,isym4,inn, iV,kj) stop 102 format(2i7,1p,e10.3) end ! ------------------------------------------------------------------ subroutine read_energies!(E,ij,isym4,inn,nen) use energylevs implicit double precision(a-h,o-y), logical(z) ! common/energylevs/ enlevs,ij,isym4,inn, vvkj, nen ! dimension E(NeN), ij(nen), inn(nen), isym4(nen),vv(nen,3),kj(nen,3) dimension kjt(3), ivt(3) ien=80 icount=0 open(unit=ien,form='formatted') do 10 i=1,nlevs !!! write(99,*)i read(ien,101,end=20) index,ijt, isymt, innt, Etemp,(iVt(j),j=1,3),& !!! 101 (KJt(j),j=1,3) !!! write (99,101) index,ijt, isymt, innt, Etemp,(iVt(j),j=1,3),& !!! !!! (KJt(j),j=1,3) !!! ij(index)=ijt isym4(index)=isymt inn(index)=innt enlevs(index)=Etemp do 30 j=1,3 kj(index,j)=kjt(j) iv(index,j)=ivt(j) 30 continue ! kj(index,1)=ij(index) icount=icount+1 10 continue 20 close(ien) write(6,*) icount," Energy levels read from stream ", ien return 101 format(i6,i4,i2,i5,f15.6,1x,3i3,2x,3i3) end ! ------------------------------------------------------------------ subroutine isym_reverse(isym,j,ipart,kmin) ! simmetrya implicit double precision(a-h,o-y), logical(z) ! if(isym .eq. 1) then if(isym .eq. 0) then ipart=0 kmin=mod(J+1,2) return endif ! if(isym .eq. 2) then if(isym .eq. 1) then ipart=0 kmin=mod(J,2) return endif ! if(isym .eq. 3) then ! ipart=1 ! kmin=mod(j,2) ! return ! endif ! if(isym .eq. 4) then ! ipart=1 ! kmin=mod(j+1,2) ! return ! endif write(6,*) "isym_reverse failed" !!! write(6,*) "isym not between 1 and 4" write(6,*) "isym not between 0 and 1" stop return end !------------------------------------------------------------------- subroutine spmain(ilev,ispe,nr,nmax1,nmax2,item,idia,gz) ! this is the effective main program of spectra implicit double precision(a-h,o-y), logical(z) common/logic/ zout, zsort, zspe, zpfun, zembed data dwl/0.0d0/,x1/1.0d0/,x0/0.0d0/ ! sort out the spectrum on frequencies if (zsort) call sortsp(nr, item, ispe) if (.not.zspe) stop rewind ispe read(ispe) zembed ! read in the nuclear spin factors read(5,101) ge, go 101 format(6f10.0) if(ge.le.x0 .and. go.le.x0) then ge= x1 go= x1 endif write(6,205) ge, go 205 format(/5x,'Even spin factor = ',f6.3,' odd spin factor = ',f6.3) ! read in the conditions for spectra required. ! read(5,101, end=92) temp, xmin, wmin, wmax, dwl, q ! print out run conditions write(6,206) temp 206 format(/5x,'Temperature set to: ',f8.2,' K'//) if (xmin.ne.x0) then write(6,207) xmin else write(6,208) endif 207 format(5x,'Minimum relative intensity required = ',f8.6) 208 format(5x,'All transitions printed out') if(wmax.ne.x0) then write(6,209) wmin, wmax, dwl else write(6,210) wmin, dwl endif 209 format(5x,'Frequency range from ',f10.3,' cm-1 to ',f10.3,' cm-1',/ & 5x,'profile half width',f10.3/) 210 format(5x,'Frequency range from ',f10.3,' cm-1 to maximum',/ & 5x,'profile half width',f10.3/) ! calculate the partition function if (zpfun) then nlev=max(nmax1,nmax2) ! call pfcalc(temp,emax,qerr,ge,go,nlev,q,ilev,idia,gz) call pfcalc(temp,ge,go,q,idia,gz) write(6,211) q 211 format(/5x,'Partition function Q =',d12.5) else if (q .le. x0) q=x1 write(6,2061) q 2061 format(/5x,' Partition function set to Q =',d12.5) endif ! calculate integrated absorption coefficients jdia= abs(idia) call spectm(xmin, wmin, wmax, dwl, & ge, go, temp, q, nr, jdia, ispe, gz) 92 continue write(6,212) 212 format(/10x,'Job completed successfully') call timer ! stop return end ! ------------------------------------------------------------------ subroutine sortsp(nr, item, ispe) ! this subroutine sorts the dipole data on ascending frequencies. ! data printed out is in cm-1, debye**2 and sec-1. ! data written to ispe for spectm is in atomic units. implicit double precision(a-h,o-y), logical(z) DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: da integer, allocatable, dimension(:,:) :: ia integer, allocatable, dimension(:) :: iperm common/logic/ zout, zsort, zspe, zpfun, zembed fmem = (64*nr)/1048576.0d0 write(6,300) fmem 300 format(/5x,'Memory required to sort transitions',f10.3,' MB') allocate(ia(nr,9), da(nr,4), iperm(nr)) if(.not.zspe) zout= .true. write(ispe) zembed ir= 0 ! read in of data from item. ! input is in atomic units. rewind item 10 read(item,end=92,err=998) ipar1, j2,kmin2,ie2,j1,kmin1,ie1, & index2,index1, e2, e1, w, s,ibase2,ibase1 ir= ir + 1 ia(ir,1)= ipar1 ia(ir,2)= j2 ia(ir,3)= kmin2 ia(ir,4)= ie2 + ibase2 - 1 ia(ir,5)= j1 ia(ir,6)= kmin1 ia(ir,7)= ie1 + ibase1 - 1 ia(ir,8)=index2 ia(ir,9)=index1 da(ir,1)= e2 da(ir,2)= e1 da(ir,3)= w da(ir,4)= s goto 10 92 continue ! now sort on frequencies ier = 0 call dpsort(da(1,3), nr, iperm, 1, ier) if(ier .ne. 0) then write(6,*) " ERROR returned by DPSORT", ier stop endif do 698 ir=1,9 call ipperm(ia(1,ir), nr, iperm, ier) if(ier .ne. 0) write(6,*) " ERROR returned by ipperm ", ier if(ir .le. 4) then call dpperm(da(1,ir), nr, iperm, ier) if(ier .ne. 0) write(6,*) " ERROR returned by dpperm ", ier endif 698 continue ! write out transitions for subroutine spectm to ispe. ! output to ispe is in atomic units. ! print out results to lineprinter if requested. ! output to lineprinter is in wavenumbers and debye. if (zout) write(6,202) 202 format(//'ipar j2 p2 i2 j1 p1 i1 index2 index1 e2 e1',9x, & 'freq s(f-i)'/) i20= 0 do 4 ir= 1,nr write(ispe) (ia(ir,ic), ic=1,9), (da(ir,ic), ic=1,4) if (zout) write(6,203) (ia(ir,ic), ic=1,9), (da(ir,ic), ic=1,4) 4 continue 203 format(i3,2x,3(1x,i3),1x,3(1x,i3),2i7,3(2x,f9.3),2(3x,e10.3)) close(unit=item) deallocate(ia, da, iperm) write(6,204) 204 format(//10x,'Sort completed successfully') call timer return 998 write(6,898) ir+1 898 format(/,5x,'Failure reading from unit item, record',i6,/) ! stop return end ! -------------------------------------------------------------- subroutine pfcalc(temp,ge,go,q,idia,gz) use energylevs implicit double precision(a-h,o-y), logical(z) ! program constants ! autocm converts atomic units to wavenumbers data x0/0.0d0/,x1/1.0d0/,& hc/ 1.9864476d-16/, & bk/ 1.3806581d-16/, & autocm/ 2.19474624d+05/ tk= bk*temp c1= hc/tk q=x0 ! calculate the partition coefficient do 3 i=1,nlevs grot=dble(2*ij(i)+1) if(abs(idia) .eq. 2) then if(isym4(i).le.2) then gg= ge*grot else gg= go*grot endif else gg=grot endif q= q + gg*exp(-(enlevs(i)-gz)*c1) 3 continue if(q.le.x0) q= x1 write(6,*) "Partition function calculated to be ", q return end ! ------------------------------------------------------------ subroutine spectm(xmin, wmin, wmax, dwl, & ge, go, temp, q, nr, jdia, ispe, gz) ! this is the main working subroutine which calculates the ! required integrated absorption coefficients. ! this routine is designed to calculate integrated absorption ! coefficients using the data produced by subroutine sortsp. ! it uses the formula: ! 4.162034*10(-19)*gg*w*(exp(-hcw"/kt) - exp(-hcw'/kt))*s(f-i) ! i(w)= ------------------------------------------------------------ ! q ! where: ! gg is the nuclear spin weighting ! times the degeneracy weighting ! w is the transition frequency ! w" is the energy of the lower level in wavenumbers ! w' is the energy of the upper level in wavenumbers ! t is the temperature at which the spectrum is required ! s(f-i) is the line strength in debye**2 for all magnetic ! components of the line. ! q is the partition function. ! this gives the intensity in cm./molecule. ! in many cases the partition function will not be supplied, and it ! is thus set at unity. this must be taken into account when using ! integrated absorption coefficients to calculate expected spectral ! intensities at given column densities. ! relative absorption coefficients are given relative to the maximum ! set to unity. ! if zemit is true, the emissivity, rather than the integrated ! absoption coefficient, is calculated from: ! (2j'+1)*gg*hcw*exp(-hcw'/kt)*aif ! j(w)= -------------------------------- ! 4pi*q ! ! New Namelist logical switches: ! Only one of the following switches may be set to true. ! ! z90 Outputs to ilist in fort.90 style, ie. index1, index2, Einstein A. ! zassign Outputs to ilist, El, frequency intensity and Einstein A along with energy level ! Assignments. ! ! use energylevs implicit double precision(a-h,o-y), logical(z) double precision, allocatable, dimension(:,:) :: a integer, allocatable, dimension(:,:) :: iqnum namelist /spe/ emin1,emax1,jmax,zplot,zemit,iplot,zfreq,zeinst, & emin2,emax2,zprof,idat,zene,tinte,zlist,ilist, & zdop,z90,zassign,prthr,npoints, xmolm common/logic/ zout, zsort, zspe, zpfun, zembed common/base/ ibase1,ibase2 ! program constants data idat/19/,iplot/20/,ilist/36/, & hc/ 1.9864476d-16/, & pi/ 3.1415927d0/, & bk/ 1.3806581d-16/, & ! autocm converts atomic units to wavenumbers autocm/ 2.19474624d+05/, & ! autode converts atomic units to debye autode/ 2.5417662d0/, & ! detosc converts from s(f-i) in debye**2 to seconds**(-1) detosc/ 3.136186d-07/, & conv / 4.162034d-19/ ! detocg converts from debye to c.g.s. units. fmem = (76*nr)/1048576.0d0 write(6,300) fmem 300 format(/5x,'Memory required to generate spectrum',f10.3,' MB'/) allocate(iqnum(nr,9), a(nr,6)) 200 format(///) ! preset namelist parameters emin1= -1.0d27 emax1= +1.0d27 emin2= -1.0d27 emax2= +1.0d27 jmax= -1 tinte= +1.0d-15 zplot= .false. zlist= .false. zprof= .false. npoints=3000 zemit= .false. zfreq= .true. zeinst= .false. zene= .false. zdop= .false. z90= .false. zassign= .false. prthr = 0.1d0 xmolm = 18.0d0 !molecular mass default (H20). ! read in namelist parameters read(5,spe) if(zeinst) zemit=.true. write(6,*) "SPECTM output requests:" if (zplot) then write(6,1010) iplot 1010 format(5x,'Data for plotting written to stream IPLOT =',i3) open(unit=iplot,form='formatted') if (zprof) then open(unit=idat,form='formatted') write(6,1020) npoints,idat 1020 format(5x,'Full line profile requested at',i5,' points,',& ' scatch IDAT =',i3) endif else write(6,1030) 1030 format(5x,'Plot file not requested') endif if (zlist) then write(6,1040) ilist 1040 format(5x,'Full line list written to stream ILIST =',i3) open(unit=ilist,form='formatted') else write(6,1050) 1050 format(5x,'Full line list not requested') endif if (zlist) write(6,1060) prthr 1060 format(/5x,'Linelist printed using relative threshold PRTHR =', & d9.2) ! determine the spectral range required nskip= 0 ntrans= 0 nread= 0 rewind ispe read(ispe) 10 read(ispe,end=90,err=99) ipar, j2, k2, ie2, j1, k1, ie1, index2,index1,e2, e1, w nread= nread+1 w= w*autocm if(w.lt.wmin) then nskip= nskip + 1 goto 10 endif if(wmax.gt.0.0) then if(w.le.wmax) then ntrans= ntrans + 1 else goto 90 endif else ntrans= ntrans + 1 endif goto 10 90 continue ! skip unrequired transitions rewind ispe read(ispe) do 1 ir=1,nskip read(ispe) 1 continue ! read in required data from ispe: do 2 ir=1,ntrans read(ispe) (iqnum(ir,ic),ic=1,9), (a(ir,ic),ic=1,4) a(ir,1)= (a(ir,1)*autocm)-gz a(ir,2)= (a(ir,2)*autocm)-gz a(ir,3)= a(ir,3)*autocm a(ir,4)= a(ir,4)*autode*autode 2 continue ! start the calculation c1= hc/(temp*bk) if (zemit) goto 40 ! calculate integrated absorption coefficients c2= conv/q xmax= 0.0d0 xint= 0.0d0 if(jdia.eq.2) then do 3 ir=1,ntrans if(iqnum(ir,1).eq.0.or.iqnum(ir,1).eq.2) then gg= ge else gg= go endif acur= c2*gg*a(ir,3)*a(ir,4)* & (exp(-a(ir,2)*c1) - exp(-a(ir,1)*c1)) a(ir,5)= acur xint= xint + acur if(acur.gt.xmax) xmax= acur 3 continue else do 4 ir=1,ntrans acur= c2*a(ir,3)*a(ir,4)* & (exp(-a(ir,2)*c1) - exp(-a(ir,1)*c1)) a(ir,5)= acur xint= xint + acur if(acur.gt.xmax) xmax= acur 4 continue endif goto 41 ! calculate emissivity coefficients 40 c2= hc*detosc/(4.0d0*pi*q) xmax= 0.0d0 xint= 0.0d0 if(jdia.eq.2) then do 43 ir=1,ntrans w= a(ir,3) if(iqnum(ir,1).eq.0.or.iqnum(ir,1).eq.2) then gg= ge else gg= go endif if(.not. zeinst) then acur= w*w*w*w*a(ir,4)*c2*gg* & exp(-a(ir,1)*c1) else acur= w*w*w*gg*a(ir,4)*detosc endif a(ir,5)= acur xint= xint + acur if(acur.gt.xmax) xmax= acur 43 continue else do 44 ir=1,ntrans w= a(ir,3) if(.not. zeinst) then acur= w*w*w*w*a(ir,4)*c2* & exp(-a(ir,1)*c1) else acur= w*w*w*a(ir,4)*detosc endif a(ir,5)= acur xint= xint + acur if(acur.gt.xmax) xmax= acur 44 continue endif 41 continue ! relative absorption or emissivity coefficients inc= 0 neg= 0 xmax1= 0.0d0 do 5 ir= 1,ntrans if(jmax.gt.-1.and.iqnum(ir,5).gt.jmax) goto 5 if(a(ir,2).lt.emin1) goto 5 if(a(ir,2).gt.emax1) goto 5 if(a(ir,1).lt.emin2) goto 5 if(a(ir,1).gt.emax2) goto 5 if(a(ir,5).gt.xmax1) xmax1= a(ir,5) 5 continue do 6 ir= 1,ntrans a(ir,6)= a(ir,5)/xmax1 if(a(ir,6).lt.xmin) goto 7 if(jmax.gt.-1.and.iqnum(ir,5).gt.jmax) goto 7 if(a(ir,2).lt.emin1) goto 7 if(a(ir,2).gt.emax1) goto 7 if(a(ir,1).lt.emin2) goto 7 if(a(ir,1).gt.emax2) goto 7 inc= inc + 1 goto 6 7 continue neg= neg + 1 6 continue if(jmax.eq.-1) then jm=500 else jm= jmax endif if(.not.zemit) then write(6,2021) xmax,xint,emin1,emax1,emin2,emax2,jm,ntrans, & inc,neg else write(6,2022) xmax,xint,emin1,emax1,emin2,emax2,jm,ntrans, & inc,neg endif 2021 format(/5x,'Maximum absorption coefficient = ',e15.6,/, & 10x,'units are cm./molecule.'//, & 5x,'Integrated band absorption =',e15.6,//, & 5x,'Transitions with lower state energy = ',e13.6,' cm-1', & ' to ',e13.6,' cm-1 considered.',/, & 5x,'transitions with upper state energy = ',e13.6,' cm-1', & ' to ',e13.6,' cm-1 considered.',/, & 5x,'Maximum allowed value of J" =',i4,//, & i11,' transitions within spectral range',/, & i11,' transitions included',/, & i11,' transitions neglected') 2022 format(/5x,'Maximum emission coefficient = ',e15.6,/, & 10x,'units are erg/molecule/sr.'//, & 5x,'Integrated band emission =',e15.6,//, & 5x,'Transitions with lower state energy = ',e13.6,' cm-1', & ' to ',e13.6,' cm-1 considered.',/, & 5x,'Transitions with upper state energy = ',e13.6,' cm-1', & ' to ',e13.6,' cm-1 considered.',/, & 5x,'maximum allowed value of j" =',i4,//, & i11,' transitions within spectral range',/, & i11,' transitions included',/, & i11,' transitions neglected') ! final output write(6,*) i20= 0 write(6,203) 203 format(/'ipar j2 p2 i2 j1 p1 i1 index2 index1 e2 e1 ', & ' freq s(f-i) abs i(w) rel i(w) a(if) '/) 204 format(i3,2x,2i3,i4,1x,2i3,i4,2i7,3f13.6,1x,1p,e16.8,3e11.3) do 8 ir=1,ntrans if(jmax.gt.-1.and.iqnum(ir,5).gt.jmax) goto 8 if(a(ir,2).lt.emin1) goto 8 if(a(ir,2).gt.emax1) goto 8 if(a(ir,1).lt.emin2) goto 8 if(a(ir,1).gt.emax2) goto 8 if(a(ir,6).lt.xmin) goto 8 w= a(ir,3) aa= w*w*w*a(ir,4)*detosc/dble(2*iqnum(ir,2) + 1) iqnum(ir,3) = abs(iqnum(ir,3) - 1) iqnum(ir,6) = abs(iqnum(ir,6) - 1) if (zlist) then if(z90) then write(ilist,208) iqnum(ir,8), iqnum(ir,9),aa else if(zassign) then ind2=iqnum(ir,8) ind1=iqnum(ir,9) write(ilist,209) isym4(ind2),iqnum(ir,4),(iv(ind2,ic),ic=1,3), & (kj(ind2,ic),ic=1,3),isym4(ind1),iqnum(ir,7),(iv(ind1,ic),ic=1,3),& (kj(ind1,ic),ic=1,3),ind2,ind1, & a(ir,2),a(ir,3),a(ir,5),aa else write(ilist,204) (iqnum(ir,ic),ic=1,9),(a(ir,ic),ic=1,6),aa endif endif endif if(a(ir,6).ge.prthr) then write(6,204) (iqnum(ir,ic),ic=1,9),(a(ir,ic),ic=1,6),aa i20= i20 + 1 endif if(zplot.and..not.zprof) then if( a(ir,6) .ge.tinte) then if ( zfreq ) then xval= a(ir,3) else xval= 10000.0/a(ir,3) endif if ( .not.zene ) then write(iplot,207) xval, a(ir,6) else write(iplot,206) xval, a(ir,6),a(ir,1),a(ir,2), & iqnum(ir,2),iqnum(ir,5) endif endif endif if(i20.eq.20) then write(6,200) i20= 0 endif 8 continue if(zplot.and.zprof) then do 9 ir=1,ntrans if( zfreq ) then a(ir,6)= a(ir,3) else a(ir,6)= 10000.0/a(ir,3) endif write(idat,*) a(ir,6), a(ir,5) 9 continue fmin= a(1,6) fmax= a(ntrans,6) call profil(idat,dwl,temp,iplot,fmin,fmax,npoints,& xmolm,zdop,zfreq) ! call profil(idat,dwl,iplot,fmin,fmax,npoints) endif 206 format(1x,f11.5,2x,f13.10,2x,f12.4,2x,f12.4,2x,i3,2x,i3) 207 format(1x,f11.5,2x,f20.10) write(6,200) write(6,205) 205 format(10x,'Spectrum completed successfully') return 99 write(6,*) nread, j2,k2,ie2,j1,k1,ie1,e2,e1,w deallocate(iqnum, a) ! stop return 208 format(2i7,1p,e10.3) 209 format(i2,i5,6i3,2x,i2,i5,6i3,2x,2i7,2x,2f11.4,1p,2e10.3) end ! ---------------------------------------------------------- subroutine getrow(row,nrow,iunit) implicit double precision (a-h,o-y) dimension row(nrow) read(iunit,err=999) row return 999 write(6,899) iunit 899 format(/,5x,'fell over in getrow reading unit = ',i6,/) stop end ! ------------------------------------------------------------ subroutine outrow(row,nrow,iunit) implicit double precision (a-h,o-y) dimension row(nrow) write(iunit) row return end subroutine timer ! prints current cpu time usage #030 implicit double precision (a-h,o-y) common/timing/itime0 write(6,10) call SYSTEM_CLOCK(itime2,irate2,imax2) itime=(itime2-itime0)/irate2 write(6,1)itime 1 format(/i10,' secs CPU time used'/) return 10 format(/) end ! -------------------------------------------------------- subroutine profil(idat,dwl,temp,iplot,fmin,fmax,npoints,& xmolm,zdop,zfreq) ! This routine calculates a Gaussian or thermal Doppler line profile ! and applies it to the previously calculated integrated absorption ! coefficient or the emissivity. The frequency dependent total ! absorption coefficient and emissivity for the molecule are ! calculated at each frequency or wavelength point. ! ! Output is to the unit idat. ! Units of the line profiled absorption coefficient are [cm**2 /molecule] ! Units of the line profiled emissivity are [ergs cm /molecule /str]. ! If units of [ergs /micron /molecule /str] are desired, the user ! must multiply by (10000.0 / (lambda[microns])**2) ! ! ARGUMENTS ! idat I/O unit, on which frequency/intensity data is stored. ! dwl Gaussian half width as specified by the user ! temp Temperature ! ipolt I/O unit, on which data is to be output. ! fmin minimum of frequency range ! fmax maximum of frequency range ! ! Input arguments from namelist SPE. ! ! npoints Number of frequency points, default is 3000. ! Zdop logical parameter, if true then a thermal Doppler ! line half width is used. The default is false, in which ! case the user specified value of the half width at half ! maximum of a Gaussian is used. ! Zfreq Default .true. units of wavenumber, and calculation ! proceeds as discussed below. If .false. units of ! wavelength in microns. For consistency of intensity units ! the line profiling is performed only in wavenumbers (cm**-1). ! But intensity data is output as a function of wavelength. ! xmolm Mean molecular mass of the species, default is 18.0d0 (H20). ! ! OTHER INPUT ! !Input from unit idat ! ! w Wavenumber (cm**-1) or wavelength (microns) of line. ! x Integrated absorption coefficient (cm / molecule) or ! emissivity (ergs /s /str) of line. ! ! OUTPUT to unit iplot ! ! wl(npoints) Wavenumber (cm**-1) or wavelength (microns) of point. ! f(npoints) Absorption cross section (cm**2 / molecule) or ! emissivity (ergs cm /s /str) at the point. ! ! The Doppler profile is: ! ! 1 ( m*c*c ) ( -m*c*c ) ! DLP = --- * sqrt(---------) * exp(----------- * Dv*Dv ) ! v ( 2*pi*k*T) ( 2*k*T*v*v ) ! ! Where: v is the central line frequency or wavenumber ! dv is the displacement from the line centre ! m is the mass of particle, T is temperature. ! ! The Doppler half width at half maximum is thus: ! ! ( 2*log(2)*k*T ) ! HW = v * sqrt( ------------ ) ! ( m*c*c ) ! ! Note: log(2) is the natural log of 2, ie FORTRAN notation. ! ! the profile can now be written: ! ! 1 ( log(2) ) ( -Dv*Dv*log(2) ) ! DLP = ---- * sqrt( -------- ) * exp( ------------- ) ! HW ( pi ) ( HW*HW ) ! ! the factor sqrt(log(2)) is often dropped from the definition of ! the Doppler half width, so simplifying the all the expressions. ! ! Here, however, we use the rigorous definition of ! half width at half maximum of the Gaussian. implicit double precision (a-h,o-y), logical(z) double precision Temp, xmolm, xln2, rtln2, rtln2pi, xlambda DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: f, wl allocate(f(npoints), wl(npoints)) write(6,*) write(6,*) "LINE PROFILING" if(zdop) then write(6,*) "Thermal Doppler line profiling requested." write(6,*) "Molecular mass set to", xmolm else write(6,*) "Gaussian profiling requested." write(6,*) "Half width at half maximum set to ", dwl, " cm**-1" endif if (zfreq) then write(6,*) "Wavenumber units of cm**-1 selected" else write(6,*) "Wavelength units of microns selected" endif df= (fmax - fmin)/dble(npoints-1) ! some constants. pi= acos(-1.0d0) xln2=log(2.0d0) rtln2=sqrt(xln2) rtln2pi=rtln2/sqrt(pi) ! zero the intensity array do 1 i= 1,npoints f(i)= 0.0d0 wl(i)= fmin + dble(i-1)*df 1 continue ! read in data from idat rewind idat 10 read(idat,*,end=11) w, x if(zfreq .eqv. .false.) then xlambda=w w=10000.0d0/w endif ! set the half width at half maximum [cm**-1] if(zdop) then hw=dop_half(temp,w,xmolm) else hw=dwl endif fac= rtln2pi/hw ! Normalisation factor [cm] do 12 i=1,npoints ! set displacement from line centre in [cm**-1] if(zfreq) then delw= wl(i) - w else delw= (10000.0d0/wl(i))-w endif if(abs(delw).gt.10.0d0*hw) goto 12 xarg= delw/hw xx= x*exp(-xarg*xarg*xln2) ! Guassian line profile, .... f(i)= f(i) + (xx*fac) ! normalise and bin. 12 continue goto 10 ! output the spectrum 11 continue do 20 i=1,npoints write(iplot,*) wl(i), f(i) 20 continue deallocate(f, wl) return end ! ---------------------------------------------------------- function dop_half(T,nu_cm, molm) ! This function calculates the Gaussian thermal Doppler half ! width at half maximum. ! INPUT: ! T: temperature ! nu_cm: Wavenumber (cm-1). ! molm: Molecular mass ! OUTPUT: ! dop_half: doppler half width (cm-1). implicit none double precision T, Nu_cm,molm, dop_half, C data C /3.5682149d-7/ dop_half=C*sqrt(T/Molm)*nu_cm return end ! ------------------------------------------------------------- function qint_H2O(T) ! This function calculates the ro-vibrational H2O partition funtion from the fit of ! Vidler and Tennyson. implicit none double precision qint_H2O, T, AZ_mol(7) double precision logZ, logtheta integer i data (AZ_mol(i),i=1,7) /-14.087469157d0,37.924324853d0, & -42.681797873d0,25.330244851d0, & -8.108512629d0,1.331068717d0, & -0.087298105d0/ logtheta=log10(T) logZ=0.0d0 do 10 i=1, 7 logZ=logZ+(az_mol(i)*(logtheta**(i-1))) 10 continue qint_H2O=10.0d0**(logZ) return end ! ------------------------------------------------------ !*DECK DPPERM SUBROUTINE DPPERM (DX, N, IPERM, IER) !C***BEGIN PROLOGUE DPPERM !C***PURPOSE Rearrange a given array according to a prescribed !C permutation vector. !C***LIBRARY SLATEC !C***CATEGORY N8 !C***TYPE DOUBLE PRECISION (SPPERM-S, DPPERM-D, IPPERM-I, HPPERM-H) !C***KEYWORDS PERMUTATION, REARRANGEMENT !C***AUTHOR McClain, M. A., (NIST) !C Rhoads, G. S., (NBS) !C***DESCRIPTION !C !C DPPERM rearranges the data vector DX according to the !C permutation IPERM: DX(I) <--- DX(IPERM(I)). IPERM could come !C from one of the sorting routines IPSORT, SPSORT, DPSORT or !C HPSORT. !C !C Description of Parameters !C DX - input/output -- double precision array of values to be !C rearranged. !C N - input -- number of values in double precision array DX. !C IPERM - input -- permutation vector. !C IER - output -- error indicator: !C = 0 if no error, !C = 1 if N is zero or negative, !C = 2 if IPERM is not a valid permutation. !C !C***REFERENCES (NONE) !C***ROUTINES CALLED XERMSG !C***REVISION HISTORY (YYMMDD) !C 901004 DATE WRITTEN !C 920507 Modified by M. McClain to revise prologue text. !C***END PROLOGUE DPPERM INTEGER N, IPERM(*), I, IER, INDX, INDX0, ISTRT DOUBLE PRECISION DX(*), DTEMP !C***FIRST EXECUTABLE STATEMENT DPPERM IER=0 IF(N.LT.1)THEN IER=1 ! CALL XERMSG ('SLATEC', 'DPPERM',& ! 'The number of values to be rearranged, N, is not positive.',IER, 1) RETURN ENDIF !C !C CHECK WHETHER IPERM IS A VALID PERMUTATION !C DO 100 I=1,N INDX=ABS(IPERM(I)) IF((INDX.GE.1).AND.(INDX.LE.N))THEN IF(IPERM(INDX).GT.0)THEN IPERM(INDX)=-IPERM(INDX) GOTO 100 ENDIF ENDIF IER=2 !!!! CALL XERMSG ('SLATEC', 'DPPERM',& !!!! 'The permutation vector, IPERM, is not valid.', IER, 1) RETURN 100 CONTINUE !C !C REARRANGE THE VALUES OF DX !C !C USE THE IPERM VECTOR AS A FLAG. !C IF IPERM(I) > 0, THEN THE I-TH VALUE IS IN CORRECT LOCATION !C DO 330 ISTRT = 1 , N IF (IPERM(ISTRT) .GT. 0) GOTO 330 INDX = ISTRT INDX0 = INDX DTEMP = DX(ISTRT) 320 CONTINUE IF (IPERM(INDX) .GE. 0) GOTO 325 DX(INDX) = DX(-IPERM(INDX)) INDX0 = INDX IPERM(INDX) = -IPERM(INDX) INDX = IPERM(INDX) GOTO 320 325 CONTINUE DX(INDX0) = DTEMP 330 CONTINUE !C RETURN END !*DECK DPSORT SUBROUTINE DPSORT (DX, N, IPERM, KFLAG, IER) !C***BEGIN PROLOGUE DPSORT !C***PURPOSE Return the permutation vector generated by sorting a given !!C array and, optionally, rearrange the elements of the array. !C The array may be sorted in increasing or decreasing order. !C A slightly modified quicksort algorithm is used. !C***LIBRARY SLATEC !C***CATEGORY N6A1B, N6A2B !C***TYPE DOUBLE PRECISION (SPSORT-S, DPSORT-D, IPSORT-I, HPSORT-H) !C***KEYWORDS NUMBER SORTING, PASSIVE SORTING, SINGLETON QUICKSORT, SORT !C***AUTHOR Jones, R. E., (SNLA) !C Rhoads, G. S., (NBS) !C Wisniewski, J. A., (SNLA) !C***DESCRIPTION !C !C DPSORT returns the permutation vector IPERM generated by sorting !C the array DX and, optionally, rearranges the values in DX. DX may !!C be sorted in increasing or decreasing order. A slightly modified !C quicksort algorithm is used. !C !C IPERM is such that DX(IPERM(I)) is the Ith value in the !C rearrangement of DX. IPERM may be applied to another array by !C calling IPPERM, SPPERM, DPPERM or HPPERM. !C !C The main difference between DPSORT and its active sorting equivalent !C DSORT is that the data are referenced indirectly rather than !C directly. Therefore, DPSORT should require approximately twice as !C long to execute as DSORT. However, DPSORT is more general. !C !C Description of Parameters !C DX - input/output -- double precision array of values to be !C sorted. If ABS(KFLAG) = 2, then the values in DX will be !C rearranged on output; otherwise, they are unchanged. !C N - input -- number of values in array DX to be sorted. !C IPERM - output -- permutation array such that IPERM(I) is the !C index of the value in the original order of the !C DX array that is in the Ith location in the sorted !C order. !C KFLAG - input -- control parameter: !C = 2 means return the permutation vector resulting from !C sorting DX in increasing order and sort DX also. !C = 1 means return the permutation vector resulting from !C sorting DX in increasing order and do not sort DX. !C = -1 means return the permutation vector resulting from !C sorting DX in decreasing order and do not sort DX. !C = -2 means return the permutation vector resulting from !C sorting DX in decreasing order and sort DX also. !C IER - output -- error indicator: !C = 0 if no error, !C = 1 if N is zero or negative, !C = 2 if KFLAG is not 2, 1, -1, or -2. !C***REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm !C for sorting with minimal storage, Communications of !C the ACM, 12, 3 (1969), pp. 185-187. !C***ROUTINES CALLED XERMSG !C***REVISION HISTORY (YYMMDD) !C 761101 DATE WRITTEN !C 761118 Modified by John A. Wisniewski to use the Singleton !C quicksort algorithm. !C 870423 Modified by Gregory S. Rhoads for passive sorting with the !C option for the rearrangement of the original data. !C 890619 Double precision version of SPSORT created by D. W. Lozier. !C 890620 Algorithm for rearranging the data vector corrected by R. !C Boisvert. !C 890622 Prologue upgraded to Version 4.0 style by D. Lozier. !C 891128 Error when KFLAG.LT.0 and N=1 corrected by R. Boisvert. !C 920507 Modified by M. McClain to revise prologue text. !C 920818 Declarations section rebuilt and code restructured to use !C IF-THEN-ELSE-ENDIF. (SMR, WRB) !C***END PROLOGUE DPSORT !C .. Scalar Arguments .. INTEGER IER, KFLAG, N !C .. Array Arguments .. DOUBLE PRECISION DX(*) INTEGER IPERM(*) !C .. Local Scalars .. DOUBLE PRECISION R, TEMP INTEGER I, IJ, INDX, INDX0, ISTRT, J, K, KK, L, LM, LMT, M, NN !C .. Local Arrays .. INTEGER IL(21), IU(21) !C .. External Subroutines .. !!!!!! EXTERNAL XERMSG !C .. Intrinsic Functions .. INTRINSIC ABS, INT !C***FIRST EXECUTABLE STATEMENT DPSORT IER = 0 NN = N IF (NN .LT. 1) THEN IER = 1 !!!!! CALL XERMSG ('SLATEC', 'DPSORT',& !!!!!! 'The number of values to be sorted, N, is not positive.', IER, 1) RETURN ENDIF !C KK = ABS(KFLAG) IF (KK.NE.1 .AND. KK.NE.2) THEN IER = 2 !!!!! CALL XERMSG ('SLATEC', 'DPSORT', & !!!! 'The sort control parameter, KFLAG, is not 2, 1, -1, or -2.', IER, 1) RETURN ENDIF !C !C Initialize permutation vector !C DO 10 I=1,NN IPERM(I) = I 10 CONTINUE !C !C Return if only one value is to be sorted !C IF (NN .EQ. 1) RETURN !C !C Alter array DX to get decreasing order if needed !C IF (KFLAG .LE. -1) THEN DO 20 I=1,NN DX(I) = -DX(I) 20 CONTINUE ENDIF !C !C Sort DX only !C M = 1 I = 1 J = NN R = .375D0 !C 30 IF (I .EQ. J) GO TO 80 IF (R .LE. 0.5898437D0) THEN R = R+3.90625D-2 ELSE R = R-0.21875D0 ENDIF !C 40 K = I !C !C Select a central element of the array and save it in location L !C IJ = I + INT((J-I)*R) LM = IPERM(IJ) !C !C If first element of array is greater than LM, interchange with LM !C IF (DX(IPERM(I)) .GT. DX(LM)) THEN IPERM(IJ) = IPERM(I) IPERM(I) = LM LM = IPERM(IJ) ENDIF L = J !C !C If last element of array is less than LM, interchange with LM !C IF (DX(IPERM(J)) .LT. DX(LM)) THEN IPERM(IJ) = IPERM(J) IPERM(J) = LM LM = IPERM(IJ) !C !C If first element of array is greater than LM, interchange !C with LM !C IF (DX(IPERM(I)) .GT. DX(LM)) THEN IPERM(IJ) = IPERM(I) IPERM(I) = LM LM = IPERM(IJ) ENDIF ENDIF GO TO 60 50 LMT = IPERM(L) IPERM(L) = IPERM(K) IPERM(K) = LMT !C !C Find an element in the second half of the array which is smaller !C than LM !C 60 L = L-1 IF (DX(IPERM(L)) .GT. DX(LM)) GO TO 60 !C !C Find an element in the first half of the array which is greater !C than LM !C 70 K = K+1 IF (DX(IPERM(K)) .LT. DX(LM)) GO TO 70 !C !C Interchange these elements !C IF (K .LE. L) GO TO 50 !C !C Save upper and lower subscripts of the array yet to be sorted !C IF (L-I .GT. J-K) THEN IL(M) = I IU(M) = L I = K M = M+1 ELSE IL(M) = K IU(M) = J J = L M = M+1 ENDIF GO TO 90 !C !C Begin again on another portion of the unsorted array !C 80 M = M-1 IF (M .EQ. 0) GO TO 120 I = IL(M) J = IU(M) !C 90 IF (J-I .GE. 1) GO TO 40 IF (I .EQ. 1) GO TO 30 I = I-1 !C 100 I = I+1 IF (I .EQ. J) GO TO 80 LM = IPERM(I+1) IF (DX(IPERM(I)) .LE. DX(LM)) GO TO 100 K = I !C 110 IPERM(K+1) = IPERM(K) K = K-1 IF (DX(LM) .LT. DX(IPERM(K))) GO TO 110 IPERM(K+1) = LM GO TO 100 !C !C Clean up !C 120 IF (KFLAG .LE. -1) THEN DO 130 I=1,NN DX(I) = -DX(I) 130 CONTINUE ENDIF !C !C Rearrange the values of DX if desired !C IF (KK .EQ. 2) THEN !C !C Use the IPERM vector as a flag. !C If IPERM(I) < 0, then the I-th value is in correct location !C DO 150 ISTRT=1,NN IF (IPERM(ISTRT) .GE. 0) THEN INDX = ISTRT INDX0 = INDX TEMP = DX(ISTRT) 140 IF (IPERM(INDX) .GT. 0) THEN DX(INDX) = DX(IPERM(INDX)) INDX0 = INDX IPERM(INDX) = -IPERM(INDX) INDX = ABS(IPERM(INDX)) GO TO 140 ENDIF DX(INDX0) = TEMP ENDIF 150 CONTINUE !C !C Revert the signs of the IPERM values !C DO 160 I=1,NN IPERM(I) = -IPERM(I) 160 CONTINUE !C ENDIF !C RETURN END !*DECK IPPERM SUBROUTINE IPPERM (IX, N, IPERM, IER) !C***BEGIN PROLOGUE IPPERM !C***PURPOSE Rearrange a given array according to a prescribed !C permutation vector. !C***LIBRARY SLATEC !C***CATEGORY N8 !C***TYPE INTEGER (SPPERM-S, DPPERM-D, IPPERM-I, HPPERM-H) !C***KEYWORDS APPLICATION OF PERMUTATION TO DATA VECTOR !C***AUTHOR McClain, M. A., (NIST) !C Rhoads, G. S., (NBS) !C***DESCRIPTION !C !C IPPERM rearranges the data vector IX according to the !C permutation IPERM: IX(I) <--- IX(IPERM(I)). IPERM could come !C from one of the sorting routines IPSORT, SPSORT, DPSORT or !C HPSORT. !C !C Description of Parameters !C IX - input/output -- integer array of values to be rearranged. !C N - input -- number of values in integer array IX. !C IPERM - input -- permutation vector. !C IER - output -- error indicator: !C = 0 if no error, !C = 1 if N is zero or negative, !C = 2 if IPERM is not a valid permutation. !C !C***REFERENCES (NONE) !C***ROUTINES CALLED XERMSG !!C***REVISION HISTORY (YYMMDD) !C 900618 DATE WRITTEN !C 920507 Modified by M. McClain to revise prologue text. !C***END PROLOGUE IPPERM INTEGER IX(*), N, IPERM(*), I, IER, INDX, INDX0, ITEMP, ISTRT !C***FIRST EXECUTABLE STATEMENT IPPERM IER=0 IF(N.LT.1)THEN IER=1 !!!!!! CALL XERMSG ('SLATEC', 'IPPERM',& !!!!!! 'The number of values to be rearranged, N, is not positive.', IER, 1) RETURN ENDIF !C !C CHECK WHETHER IPERM IS A VALID PERMUTATION !C DO 100 I=1,N INDX=ABS(IPERM(I)) IF((INDX.GE.1).AND.(INDX.LE.N))THEN IF(IPERM(INDX).GT.0)THEN IPERM(INDX)=-IPERM(INDX) GOTO 100 ENDIF ENDIF IER=2 !!!!!! CALL XERMSG ('SLATEC', 'IPPERM',& !!!!!! 'The permutation vector, IPERM, is not valid.', IER, 1) RETURN 100 CONTINUE !C !C REARRANGE THE VALUES OF IX !C !C USE THE IPERM VECTOR AS A FLAG. !C IF IPERM(I) > 0, THEN THE I-TH VALUE IS IN CORRECT LOCATION !C DO 330 ISTRT = 1 , N IF (IPERM(ISTRT) .GT. 0) GOTO 330 INDX = ISTRT INDX0 = INDX ITEMP = IX(ISTRT) 320 CONTINUE IF (IPERM(INDX) .GE. 0) GOTO 325 IX(INDX) = IX(-IPERM(INDX)) INDX0 = INDX IPERM(INDX) = -IPERM(INDX) INDX = IPERM(INDX) GOTO 320 325 CONTINUE IX(INDX0) = ITEMP 330 CONTINUE !C RETURN END