diff --git a/src/radex_src/commondata.f90 b/src/radex_src/commondata.f90 index 555314b..7dfc279 100644 --- a/src/radex_src/commondata.f90 +++ b/src/radex_src/commondata.f90 @@ -3,7 +3,7 @@ MODULE CommonData IMPLICIT NONE !file for input and output - character(200) :: outfile,molfile,specref + character(200) :: outfile,molfile,molfileold,specref character(*), PARAMETER :: radat='./data/' character(*), PARAMETER :: version = '30nov2011' character(*), PARAMETER :: logfile = './radex.log' @@ -115,8 +115,12 @@ MODULE CommonData ! ctot: total collision rate ! xpop: level populations - + INTEGER :: id(maxpart) ! to identify collision partners + INTEGER :: lcu(maxpart,maxcoll),lcl(maxpart,maxcoll) + REAL(dp) :: coll(maxpart,maxcoll,maxtemp) + REAL(dp) :: temp(maxpart,maxtemp) ! collision temperatures + character*120 collref ! text about source of collisional data !For development / maintenance purposes: LOGICAL, PARAMETER :: debug =.False. -END MODULE CommonData \ No newline at end of file +END MODULE CommonData diff --git a/src/radex_src/io.f90 b/src/radex_src/io.f90 index 6e41636..3231932 100644 --- a/src/radex_src/io.f90 +++ b/src/radex_src/io.f90 @@ -16,19 +16,7 @@ SUBROUTINE ReadData(success) INTEGER :: icoll ! to loop over collisional transitions INTEGER :: dummy ! to skip part of the file - INTEGER :: id(maxpart) ! to identify collision partners - - ! upper/lower levels of collisional transition - INTEGER :: lcu(maxcoll),lcl(maxcoll) - REAL(dp) :: coll(maxpart,maxcoll,maxtemp) - REAL(dp) :: colld(maxpart,maxlev,maxlev) - REAL(dp) :: ediff - REAL(dp) :: temp(maxtemp) ! collision temperatures - - character*120 collref ! text about source of collisional data - ! to interpolate rate coeffs - INTEGER :: iup,ilo,nint REAL(dp) :: tupp,tlow,fint ! to verify matching of densities and rates @@ -38,6 +26,7 @@ SUBROUTINE ReadData(success) REAL(dp) :: opr ! Executable part begins here. + OPEN(unit=11,file=molfile,status='old',err=99) ! in the header, every second line is a comment 101 format(a) @@ -107,44 +96,79 @@ SUBROUTINE ReadData(success) IF (ntemp.lt.0) stop 'error: too few collision temperatures defined' IF (ntemp.gt.maxtemp) stop 'error: too many collision temperatures' READ(11,*) - READ(11,*) (temp(itemp),itemp=1,ntemp) + READ(11,*) (temp(ipart,itemp),itemp=1,ntemp) READ(11,*) IF (debug) write(*,*) 'ready to read ',ncoll,' rates for '& &,ntemp,' temperatures for partner ',ipart DO icoll=1,ncoll - READ(11,*) dummy,lcu(icoll),lcl(icoll),& + READ(11,*) dummy,lcu(ipart,icoll),lcl(ipart,icoll),& &(coll(ipart,icoll,itemp),itemp=1,ntemp) IF ((dummy.lt.1).or.(dummy.gt.ncoll)) stop 'error:illegal collision number' END DO + END DO + + close(11) + RETURN + ! IF (debug) WRITE(*,*)'ctot=',(ctot(ilev),ilev=1,nlev) + 99 write(*,*) 'error opening data file ',molfile + STOP + END SUBROUTINE ReadData + + + SUBROUTINE InterpColli(success) + INTEGER, INTENT(INOUT) :: success + + ! Reads molecular data files (2003 format) + + INTEGER :: ilev,jlev ! to loop over energy levels + INTEGER :: iline ! to loop over lines + INTEGER :: ipart,jpart ! to loop over collision partners + INTEGER :: itemp ! to loop over collision temperatures + INTEGER :: icoll ! to loop over collisional transitions + INTEGER :: dummy ! to skip part of the file + + ! upper/lower levels of collisional transition + REAL(dp) :: ediff + REAL(dp) :: colld(maxpart,maxlev,maxlev) + + ! to interpolate rate coeffs + INTEGER :: iup,ilo,nint + REAL(dp) :: tupp,tlow,fint + + ! to verify matching of densities and rates + logical found + + ! to calculate thermal o/p ratio for H2 + REAL(dp) :: opr ! interpolate array coll(ncol,ntemp) to desired temperature ! Must DO this now because generally, rates with different partners ! are calculated for a different set of temperatures - + DO ipart=1,npart IF (ntemp.le.1) THEN DO icoll=1,ncoll - iup=lcu(icoll) - ilo=lcl(icoll) + iup=lcu(ipart,icoll) + ilo=lcl(ipart,icoll) colld(ipart,iup,ilo)=coll(ipart,icoll,1) END DO else - IF (tkin.gt.temp(1)) THEN - IF (tkin.lt.temp(ntemp)) THEN + IF (tkin.gt.temp(ipart,1)) THEN + IF (tkin.lt.temp(ipart,ntemp)) THEN !=== interpolation : DO itemp=1,(ntemp-1) - IF (tkin.gt.temp(itemp).and.tkin.le.temp(itemp+1)) nint=itemp + IF (tkin.gt.temp(ipart,itemp).and.tkin.le.temp(ipart,itemp+1)) nint=itemp END DO - tupp=temp(nint+1) - tlow=temp(nint) + tupp=temp(ipart,nint+1) + tlow=temp(ipart,nint) fint=(tkin-tlow)/(tupp-tlow) !cc db !c write(*,*) 'ipart,nint,fint: ',ipart,nint,fint DO icoll=1,ncoll - iup=lcu(icoll) - ilo=lcl(icoll) + iup=lcu(ipart,icoll) + ilo=lcl(ipart,icoll) colld(ipart,iup,ilo)=coll(ipart,icoll,nint)& &+fint*(coll(ipart,icoll,nint+1)& &-coll(ipart,icoll,nint)) @@ -153,25 +177,25 @@ SUBROUTINE ReadData(success) END DO else !=== Tkin too high : - ! IF (tkin.ne.temp(ntemp)) THEN + ! IF (tkin.ne.temp(ipart,ntemp)) THEN ! WRITE(*,*)' Warning : Tkin higher than temperatures '/& ! &/'for which rates are present.' ! END IF DO icoll=1,ncoll - iup=lcu(icoll) - ilo=lcl(icoll) + iup=lcu(ipart,icoll) + ilo=lcl(ipart,icoll) colld(ipart,iup,ilo)=coll(ipart,icoll,ntemp) END DO END IF else ! Tkin too low : - IF (tkin.ne.temp(1)) THEN + IF (tkin.ne.temp(ipart,1)) THEN WRITE(*,*)' Warning : Tkin lower than temperatures '//& &'for which rates are present.' END IF DO icoll=1,ncoll - iup=lcu(icoll) - ilo=lcl(icoll) + iup=lcu(ipart,icoll) + ilo=lcl(ipart,icoll) colld(ipart,iup,ilo)=coll(ipart,icoll,1) END DO END IF @@ -183,8 +207,6 @@ SUBROUTINE ReadData(success) IF (debug) write(*,*) 'readdata: rate coeffs' - close(11) - !$$$ IF (debug) THEN !$$$ WRITE(*,*)colld(1,1,1),colld(1,1,2),colld(1,1,3) !$$$ WRITE(*,*)colld(1,2,1),colld(1,2,2),colld(1,2,3) @@ -282,11 +304,7 @@ SUBROUTINE ReadData(success) !$$$ WRITE(*,*)crate(2,1),crate(2,2),crate(2,3) !$$$ WRITE(*,*)crate(3,1),crate(3,2),crate(3,3) !$$$ END IF - RETURN - ! IF (debug) WRITE(*,*)'ctot=',(ctot(ilev),ilev=1,nlev) - 99 write(*,*) 'error opening data file ',molfile - STOP - END SUBROUTINE ReadData + END SUBROUTINE InterpColli SUBROUTINE GetInputs INTEGER :: ipart ! loop over collision partners @@ -624,4 +642,4 @@ SUBROUTINE ParseInputDictionary(dictionary,success) fmax = dmin1(3.d7,fmax) END SUBROUTINE ParseInputDictionary -END MODULE IO \ No newline at end of file +END MODULE IO diff --git a/src/radex_src/solver.f90 b/src/radex_src/solver.f90 index af591ee..a0c4ce9 100644 --- a/src/radex_src/solver.f90 +++ b/src/radex_src/solver.f90 @@ -32,6 +32,7 @@ SUBROUTINE Matrix(niter,conv,success_flag) REAL(dp) :: sumx ! summed radiative rate REAL(dp) :: total ! to normalize populations REAL(dp) :: tsum,thistex ! to check convergence + REAL(dp) :: lr ! learning rate to solve LOGICAL :: conv ! are we converged? @@ -297,8 +298,13 @@ SUBROUTINE Matrix(niter,conv,success_flag) END IF !sb301111 now DO the underrelaxation! + IF(niter.le.20) THEN + lr=0.1+0.01*niter + else + lr=5.0/niter+0.05 + END IF DO ilev=1,nlev - xpop(ilev)=0.3*xpop(ilev)+0.7*xpopold(ilev) + xpop(ilev)=lr*xpop(ilev)+(1.0-lr)*xpopold(ilev) end do return @@ -311,17 +317,19 @@ FUNCTION EscProb(TAU) REAL(dp) :: EscProb,beta,tau REAL(dp) :: taur !optical radius - taur = tau/2.0 SELECT CASE (method) CASE(1) !Uniform sphere formula from Osterbrock (Astrophysics of !Gaseous Nebulae and Active Galactic Nuclei) Appendix 2 !with power law approximations for large and small tau - IF (abs(taur).lt.0.1) THEN + taur=tau/2.0 + IF (taur.lt.0) THEN + beta=1.0-tau + ELSE IF (taur.lt.0.1) THEN beta = 1.d0-0.75d0*taur+(taur**2.)/2.5d0& &-(taur**3.)/6.d0+(taur**4.)/17.5d0 - ELSE IF(abs(taur).gt.5.d1) THEN + ELSE IF(taur.gt.5.d1) THEN beta = 0.75d0/taur ELSE beta = 0.75d0/taur*(1.d0-1.d0/(2.d0*(taur**2.))+& @@ -333,25 +341,30 @@ FUNCTION EscProb(TAU) !Expanding sphere = Large Velocity Gradient (LVG) or Sobolev case. !Formula from De Jong, Boland and Dalgarno (1980, A&A 91, 68) !corrected by factor 2 in order to match EscProb(TAU=0)=1 - IF (abs(taur).lt.0.01) THEN - beta = 1.0 - else if(abs(taur).lt.7.0) THEN - beta = 2.0*(1.0 - dexp(-2.34*taur))/(4.68*taur) - else - beta = 2.0/(taur*4.0*(sqrt(log(taur/sqrt(pi))))) - END IF - + taur=2.34*tau + IF (taur.lt.0) THEN + beta=1.0-tau + ELSE IF (taur.lt.0.01) THEN + beta = 1.0-taur/2.0+(taur**2.)/6.d0 + ELSE IF(taur.gt.5.d1) THEN + beta = 1.0/taur + ELSE + beta = (1.0-dexp(-taur))/taur + ENDIF CASE (3) !Slab geometry (e.g., shocks): de Jong, Dalgarno & Chu 1975, !ApJ 199, 69 (again with power law approximations) - IF (abs(3.0*tau).lt.0.1) THEN - beta = 1.0 - 1.5*(tau + tau**2.) - ELSE IF (abs(3.0*tau).gt.50.0) THEN - beta = 1.0d0/(3.0*tau) + taur=3.0*tau + IF (taur.lt.0) THEN + beta=1.0-tau + ELSE IF (taur.lt.0.01) THEN + beta = 1.0-taur/2.0+(taur**2.)/6.d0 + ELSE IF(taur.gt.5.d1) THEN + beta = 1.0/taur ELSE - beta = (1.0d0 - dexp(-3.0*tau))/(3.0*tau) - END IF + beta = (1.0-dexp(-taur))/taur + ENDIF CASE DEFAULT WRITE(*,*) 'Error: Escape probability method undefined' STOP diff --git a/src/radex_src/wrap.f90 b/src/radex_src/wrap.f90 index a1f13bb..bf67553 100644 --- a/src/radex_src/wrap.f90 +++ b/src/radex_src/wrap.f90 @@ -44,7 +44,10 @@ SUBROUTINE from_params(molfileIn,tkinIn,tbgIn,cdmolIn,densityIn,& IF (success_flag .ne. 1) RETURN IF (DEBUG) write(*,*) 'calling readdata' - CALL ReadData(success_flag) + IF (molfile .ne. molfileold) THEN + CALL ReadData(success_flag) + END IF + CALL InterpColli(success_flag) IF (success_flag .ne. 1) RETURN ! Calculate background radiation field @@ -88,6 +91,7 @@ SUBROUTINE from_params(molfileIn,tkinIn,tbgIn,cdmolIn,densityIn,& niter=niter+1 END IF END DO + molfileold=molfile END SUBROUTINE from_params SUBROUTINE from_dict(inputDictionary,success_flag,nlines,Qup,Qlow,lineOutputs) @@ -111,7 +115,10 @@ SUBROUTINE from_dict(inputDictionary,success_flag,nlines,Qup,Qlow,lineOutputs) ! Read data file IF (success_flag .ne. 1) RETURN IF (DEBUG) write(*,*) 'calling readdata' - CALL ReadData(success_flag) + IF (molfile .ne. molfileold) THEN + CALL ReadData(success_flag) + END IF + CALL InterpColli(success_flag) IF (success_flag .ne. 1) RETURN ! Calculate background radiation field @@ -155,4 +162,5 @@ SUBROUTINE from_dict(inputDictionary,success_flag,nlines,Qup,Qlow,lineOutputs) niter=niter+1 END IF END DO + molfileold=molfile END SUBROUTINE from_dict