upload radex origin revision

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
+END MODULE CommonData

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
+END MODULE IO

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

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