program BYTe_CH4 implicit none ! integer, parameter :: ik = selected_int_kind(8) ! "Normal" integers. This must map on ! C "int" type, or the DX interface won't ! work correctly. integer, parameter :: hik = selected_int_kind(16) ! "Pointer" integers - sufficient to store ! memory address integer, parameter :: drk = selected_real_kind(12,25) ! "Double" reals and complex (complexi? :-) integer, parameter :: rk = selected_real_kind(12,25) ! "Normal" reals and complex (complexi? :-) integer, parameter :: ark = selected_real_kind(25,32) ! "Accurate" reals and complex (complexi? :-) ! integer(ik) :: isym,nsym,npoints,ipoint,j,jf,ji,ilevel,ilevelf,ileveli,igamma,igammaf,igammai,info,j0,i,ilevelf_,ileveli_ integer(ik) :: quantaf(0:10),quantai(0:10),imode,nmodes,localf(0:10),locali(0:10) !,normalf(0:10),normali(0:10) integer(ik) :: ib,ie,nfiles,j_t,v(9),ilog,ienerlow integer(ik) :: verbose=2,smax=200,nchar integer(hik):: Nintens(-60:60) ! integer(ik) :: nlevels,nn,itemp,ipartf,NNN real(rk) :: int_vs_enerlow(0:12000) real(rk) :: temp,temp0,freql,freqr,gns(10),halfwidth,meanmass,partfunc,dfreq,dfreq_,acoef,abscoef,& energy,energyf,energyi,tranfreq,cmcoef,beta,ln2,pi,dpwcoef,absthresh,emcoef,beta0,intband,tranfreq0 real(rk) :: de,x0,xp,xm,elow=100000. real(rk) :: thresh = 1e-16 real(rk),allocatable :: freq(:),intens(:),energies(:),pf(:,:) integer(ik),allocatable :: sym(:),n(:,:),jrot(:),krot(:),taurot(:),l(:,:),symv(:),symr(:),gtot(:) character(4) symlab(10),a_fmt character(1) invchar character(60) intfilename(200),enrfilename,swapfile character(30) proftype,specttype character(100),allocatable :: state_string(:) character(300) :: string_tmp character(len=200) :: my_fmt !text variable containing formats for reads/writes character(len=9) :: npoints_fmt !text variable containing formats for reads/writes ! logical :: swap = .false.,partfunc_do=.false.,conv_units = .false. real(rk),parameter :: planck=6.6260693d-27,avogno=6.0221415d+23,vellgt=2.99792458d+10,boltz=1.380658d-16,c2=1.4387751601679204d0 !read number of intensity files read*,nfiles if (nfiles>200) then print('("Too many files (>200):",i8)'),nfiles stop 'Too many files' endif !read intensities filenames do i = 1,nfiles read*,intfilename(i) !print('(1x,50a)'),intfilename(i) enddo !read energies filename read*,enrfilename !print('(1x,50a)'),enrfilename !read swapfile filename read*,swapfile !print('(1x,50a)'),swapfile !read number of characters read*,nchar write(a_fmt,'(I3)') nchar write(a_fmt,'("a",a3)') adjustl(a_fmt) ! write(my_fmt,'(a,a4,a,a4,a)') '(1x,2es16.8,1x,i4,1x,f12.4," <- ",i4,1x,f12.4,5x,a,"<- ",a))' !print*,my_fmt ! !read nuclear spin statistics weights ! read*,(symlab(isym),gns(isym),isym=1,nsym) !read initial state quantum numbers ! read*,(quantai(imode),imode=0,nmodes) !read final state quantum numbers ! read*,(quantaf(imode),imode=0,nmodes) !read temperature, partition function read*,temp,partfunc !read frequency range and number of points read*,freql,freqr,npoints write(npoints_fmt,'(i9)') npoints !read type of profiling read*,proftype !read type of spectra read*,specttype halfwidth = 0 ; meanmass = 0.0 ; absthresh = 0.0 ; ! !read gaussian half-width or molecular mean mass select case (trim(proftype)) case ('gauss','gauss-um'); read*,halfwidth case ('doppl','doppl-um'); read*,meanmass case ('stick','HITRAN','hitran','trunck'); read*,absthresh case ('rect ','bin '); read*,halfwidth halfwidth = min((freqr-freql)/real(npoints,8),halfwidth) case default ; stop 'illegal profile-type' end select !read the threshold read*,thresh !read the threshold ! read*,elow !compute constants ln2=log(2.0_rk) pi=acos(-1.0_rk) !beta=planck*vellgt/(boltz*temp) ! beta=c2/temp ! !cmcoef=avogno/(8.0d0*pi*vellgt) !emcoef=avogno*planck*vellgt/(4.0d0*pi) ! cmcoef=1.0_rk/(8.0d0*pi*vellgt) emcoef=1.0_rk*planck*vellgt/(4.0d0*pi) dfreq=(freqr-freql)/real(npoints-1,rk) ! ! half width for Doppler profiling dpwcoef=sqrt(2.0*ln2*boltz*avogno)/vellgt dpwcoef = dpwcoef*sqrt(temp/meanmass) if (proftype(6:8)=='-um') conv_units = .true. ! if (proftype(1:5)=='doppl') then if (verbose>=2) print('("alpha = ",f18.8)'),dpwcoef endif ! if (proftype(1:5)=='gauss') then ! if (dfreq>halfwidth*100.0) then ! if (verbose>=2) print('(" dnu >> alpha")') ! elseif(halfwidth>dfreq*100.0) then ! if (verbose>=2) print('(" dnu << alpha")') ! else ! if (verbose>=2) print('(" dnu ~ alpha")') ! endif endif ! !set up freq. grid allocate(freq(npoints),stat=info); if (info/=0) stop 'error: freq is out of memory' allocate(intens(npoints),stat=info); if (info/=0) stop 'error: intens is out of memory' forall(ipoint=1:npoints) freq(ipoint)=freql+real(ipoint-1,rk)*dfreq !count number of lines (levels) in the Energy files open(unit=1,file=trim(enrfilename)) i = 0 do ! read(1,"(i8)",end=10) nn ! i = i + 1 ! cycle ! 10 exit end do ! rewind(1) ! nlevels = i ! allocate(sym(nlevels),energies(nlevels),Jrot(nlevels),state_string(nlevels),gtot(nlevels),stat=info) if (info/=0) stop 'error: sym,energies,Jrot is out of memory' ! ! In case of specttype = partfunc the partition function will be computed for a series of temperatures. ! npoints stands for the number of temperature steps, and the frequency limits as temperature limits ! allocate(pf(0:3,npoints),stat=info) if (trim(specttype)=='partfunc') then ! !print*,halfwidth,thresh ! ipartf = int(halfwidth,4) ! dfreq=temp/real(npoints,rk) ! if (ipartf<0.or.ipartf>3) stop "illegal partfunc component, can be only 0,1,2,3" ! if (verbose>=1) print('("!",i1,a4,( 1x,'//npoints_fmt//'( 5x,"T=",f8.2,5x) ) )'),ipartf,' J ',(i*dfreq,i=1,npoints) ! if (info/=0) stop 'error: pf' endif ! ! prepare computing the partition function if (partfunc<=0) then partfunc=0 ; j0=0 partfunc_do = .true. pf = 0 endif ! ! start reading the energies ! do i=1,nlevels ! !read(1,*) nn,energies(i),gtot(i),jrot(i),invchar,sym(i),NNN,n(i,1:7),j_t,krot(i),taurot(i),v(1:6),symv(i) ! !read(1,*) nn,energies(i),gtot(i),jrot(i),string_tmp,state_string(i)(1:nchar) ! read(1,"(a250)") string_tmp ! read(string_tmp,*) nn,energies(i),gtot(i),jrot(i) read(string_tmp(42:),"("//a_fmt//")") state_string(i)(1:nchar) ! !print("("//a_fmt//")"),adjustl(string_tmp(42:min(nchar,250))) ! j = jrot(i) igamma = sym(i) energy = energies(i) ! if (partfunc_do) then if (j/=j0) then ! if (trim(specttype)=='partfunc') then if (verbose>=1.and.npoints<1000) print('("!",i4,3(1x,'//npoints_fmt//'es20.8))'),j0,pf(ipartf,1:npoints) else if (verbose>=2) print('("#",i4,1x,es16.8)'),j0,partfunc endif ! endif ! ! symmetry number ! partition function ! !partfunc=partfunc+real((2*j+1),rk)*gns(igamma)*exp(-beta*energy) ! partfunc=partfunc+gtot(i)*exp(-beta*energy) ! j0 = j endif ! if (trim(specttype)=='partfunc') then ! do itemp = 1,npoints ! if (energy>freqr) cycle ! temp0 = real(itemp,rk)*dfreq ! beta0 = planck*vellgt/(boltz*temp0) ! ! apart from the partition function, the first two derivatieves are also computed: ! !pf(0,itemp) = pf(0,itemp) + real((2*j+1),rk)* gns(igamma)*exp(-beta0*energy) !pf(1,itemp) = pf(1,itemp) + real((2*j+1),rk)*(beta0*energy) *gns(igamma)*exp(-beta0*energy) !pf(2,itemp) = pf(2,itemp) + real((2*j+1),rk)*(beta0*energy)**2*gns(igamma)*exp(-beta0*energy) ! pf(0,itemp) = pf(0,itemp) + gtot(i)*exp(-beta0*energy) pf(1,itemp) = pf(1,itemp) + gtot(i)*exp(-beta0*energy)*(beta0*energy) pf(2,itemp) = pf(2,itemp) + gtot(i)*exp(-beta0*energy)*(beta0*energy)**2 ! !qp = temp0*pf(1,itemp) !qpp= temp0**2*pf(2,itemp)+2.0d0*qp ! pf(3,itemp) = ( pf(2,itemp)/pf(0,itemp)-(pf(1,itemp)/pf(0,itemp))**2 ) ! enddo ! endif ! end do close(1) ! ! print out the computed part. function objects and finish ! if (trim(specttype)=='partfunc') then ! if (verbose>=1) print('("!",i4,3(1x,'//npoints_fmt//'es20.8/))'),j0,pf(ipartf,1:npoints) ! if (verbose>=1) print('("!0",i4,1x,'//npoints_fmt//'es20.8)'),j0,pf(0,1:npoints) if (verbose>=1) print('("!1",i4,1x,'//npoints_fmt//'es20.8)'),j0,pf(1,1:npoints) if (verbose>=1) print('("!2",i4,1x,'//npoints_fmt//'es20.8)'),j0,pf(2,1:npoints) if (verbose>=1) print('("!3",i4,1x,'//npoints_fmt//'es20.8)'),j0,pf(3,1:npoints) ! if (verbose>=2) print('(1x,a,1x,es16.8)'),'! partition function value is',partfunc ! if (verbose>=1) then print('(a)'),"Partition function" do itemp = 1,npoints temp0 = real(itemp,rk)*dfreq print('(f12.2,1x,es20.8)'),temp0,pf(0,itemp) enddo endif ! stop ! else ! if (verbose>=2) print('("#",i4,1x,es16.8)'),j0,partfunc ! endif ! if (verbose>=2) print('(1x,a,1x,es16.8)'),'# partition function value is',partfunc ! intens=0.0 ! ! the selected by the frequency limits transitions will be stored to the temporary file 'swapfile', ! which can be later used in stead of the Transition files - this will make the simualtions faster ! if (all(trim(swapfile)/=(/"null","NULL","none","NONE"/))) then open(unit=11,file=trim(swapfile),action='write',status='replace') swap = .true. endif ! !start loop over all transitions ! intens=0.0 intband = 0.0 ! !ji_ = 0 ; jf_ = 0 ! ilevelf_ = 0 ileveli_ = 0 ! Nintens = 0 int_vs_enerlow = 0 ! do i = 1,nfiles ! open(unit=1,file=trim(intfilename(i))) do ! read new line from intensities file ! read(1,*,end=20) ilevelf,ileveli,acoef ! !read(1,*,end=20) ileveli,ilevelf,acoef ! energyf = energies(ilevelf) energyi = energies(ileveli) ! if (energyf-energyi<0) then write(6,"(2i12,2x,3es16.8)"),ilevelf,ileveli,acoef,energyf,energyi stop 'wrong order of indices' endif ! if (elowfreqr) cycle ! ! if transition frequency is out of selected range if (tranfreq0>freqr+10.0*halfwidth.or.tranfreq0absthresh) print(my_fmt), & tranfreq,abscoef,jrot(ilevelf),energyf, jrot(ileveli),energyi,state_string(ilevelf)(1:nchar),& state_string(ileveli)(1:nchar) !if (abscoef>absthresh) print('((1x,2es16.8,1x,i4,1x,f12.4," <- ",i4,1x,f12.4,2x,my_fmt,"<=",my_fmt))'), & ! tranfreq,abscoef,jrot(ilevelf),energyf, jrot(ileveli),energyi,state_string(ilevelf)(1:nchar),state_string(ileveli)(1:nchar) !if (abscoef>absthresh) print('((1x,2es16.8,2x,a3,1x,i4,2x,4i4,1x,2i4,1x,i3,1x,a3,3x,2i4,1x,f12.4," <- ",a3,1x,i4,2x,4i4,1x,2i4,1x,i3,1x,a3,3x,2i4,1x,f12.4))'), & ! tranfreq,abscoef,symlab(sym(ilevelf)),jrot(ilevelf),n(ilevelf,1:7),symlab(symv(ilevelf)),krot(ilevelf),taurot(ilevelf),energyf,& ! symlab(sym(ileveli)),jrot(ileveli),n(ileveli,1:7),symlab(symv(ileveli)),krot(ileveli),taurot(ileveli),energyi !print('((1x,2es16.8,2x,a3,1x,i4,2x,6i4,3x,2i4,1x,f12.4," <- ",a3,1x,i4,2x,6i4,3x,2i4,1x,f12.4))'), & ! tranfreq,abscoef,symlab(sym(ilevelf)),jrot(ilevelf),n(ilevelf,1:6),krot(ilevelf),taurot(ilevelf),energyf,& ! symlab(sym(ileveli)),jrot(ileveli),n(ileveli,1:6),krot(ileveli),taurot(ileveli),energyi cycle end if ! if (abscoefabsthresh) write(1000+int(tranfreq/100),"(2i12,2x,es16.8,2x,f16.6)"),ilevelf,ileveli,acoef,tranfreq cycle end if ! ib = max(nint( ( tranfreq0-halfwidth*10.0-freql)/dfreq )+1,1) ie = min(nint( ( tranfreq0+halfwidth*10.0-freql)/dfreq )+1,npoints) ! ! normalize absorption coefficient ! !read gaussian half-width or molecular mean mass select case (trim(proftype(1:5))) ! case ('gauss','doppl') ! !abscoef=abscoef*sqrt(ln2/pi)/halfwidth ! !$omp parallel do private(ipoint,dfreq_,xp,xm,de) shared(intens) schedule(dynamic) do ipoint=ib,ie ! dfreq_=freq(ipoint)-tranfreq0 !if (abs(dfreq_)>halfwidth*10.0) cycle ! xp = sqrt(ln2)/halfwidth*(dfreq_)+x0 xm = sqrt(ln2)/halfwidth*(dfreq_)-x0 ! de = erf(xp)-erf(xm) ! intens(ipoint)=intens(ipoint)+abscoef*0.5_rk/dfreq*de ! enddo !$omp end parallel do ! !print('(2(1x,es16.8),i4,i4,f11.4,f11.4,f11.4,es16.8,i6,i6)'),freq(1),intens(1),ji,jf,energyi,energyf,tranfreq,linestr,ileveli,ilevelf ! !ji_ = ji ; jf_ = jf ! case ('rect'); ! !$omp parallel do private(ipoint,dfreq_) shared(intens) schedule(dynamic) do ipoint=ib,ie dfreq_=freq(ipoint)-tranfreq0 !if (abs(dfreq_)>halfwidth*10.0) cycle intens(ipoint)=max(intens(ipoint),abscoef) enddo !$omp end parallel do ! !print('(2(1x,es16.8),i4,i4,f11.4,f11.4,f11.4,es16.8,i6,i6)'),freq(1),intens(1),ji,jf,energyi,energyf,tranfreq,linestr,ileveli,ilevelf ! case ('bin'); ! !$omp parallel do private(ipoint,dfreq_) shared(intens) schedule(dynamic) do ipoint=ib,ie dfreq_=freq(ipoint)-tranfreq0 !if (abs(dfreq_)>halfwidth*10.0) cycle ! intens(ipoint)=intens(ipoint)+abscoef/dfreq ! enddo !$omp end parallel do ! !print('(2(1x,es16.8),i4,i4,f11.4,f11.4,f11.4,es16.8,i6,i6)'),freq(1),intens(1),ji,jf,energyi,energyf,tranfreq,linestr,ileveli,ilevelf ! end select ! cycle 20 exit enddo ! enddo close(1) if (proftype(1:5)=='doppl'.or.proftype(1:5)=='gauss'.or.proftype(1:4)=='rect'.or.proftype(1:3)=='bin') then ! !print out the generated intensities convoluted with a given profile ! print('(2(1x,es16.8))'),(freq(ipoint),intens(ipoint),ipoint=1,npoints) ! print('("Total intensity (sum):",es16.8," (int):",es16.8)'), intband,sum(intens)*dfreq ! else ! print('("Total intensity = ",es16.8)'), intband ! endif ! if (verbose>=3) print('(/"Intensity distribution: ")') ! do i = lbound(Nintens,dim=1),ubound(Nintens,dim=1) ! if (verbose>=3) print('(i4,2x,i10)'), i,Nintens(i) ! enddo ! !do i = lbound(int_vs_enerlow,dim=1),ubound(int_vs_enerlow,dim=1) ! ! ! if (verbose>=3) print('(f8.1,2x,es16.8)'), real(i),int_vs_enerlow(i) ! ! !enddo ! end program BYTe_CH4