c kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum c Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet c c This file must be used under the terms of the CeCILL license. c This source file is licensed as described in the file COPYING, which c you should have received as part of this distribution. The terms c are also available at c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt c subroutine calc_k(nu,Tref,T,index,P,Ps,nu0,mol_mass, & Sref,gamma_air,gamma_self,Elow, & n_air,delta_air,n_self, & c2,density,drho_dP,drho_dT,drho_dx, & Q,Qref,dQ_dT, & line_profile,usl,sl_choice,slwr,slam,slas, & sens_T,sens_P,sens_x, & sigma,k,dsigma_dP,dsigma_dT,dsigma_dx,dk_dP,dk_dT,dk_dx) implicit none include 'max.inc' c c Purpose: to compute the absorption coefficient c for a single spectral line, at a given wavenumber c c Inputs: c + nu: wavenumber [cm¯¹] c + Tref: reference temperature [K] c + T: temperature [K] c + index of the molecule c + P: pressure [atm] c + Ps: partial pressure of the given molecular species [atm] c + nu0: vacuum wavenumber [cm¯¹] c + mol_mass: molar mass of the molecular species [g/mol] c + Sref: line intensity [cm¯¹/molecule.cm¯²] at standard Tref c + gamma_air: air-broadened half-width [cm¯¹.atm¯¹] c + gamma_self: self-broadened half-width [cm¯¹.atm¯¹] c + Elow: lower-state energy [cm¯¹] c + n_air: temperature-dependance exponent for gamma_air [unitless] c + delta_air: air pressure-induced line shift [cm¯¹.atm¯¹] c + n_self: temperature-dependance exponent for gamma_self [unitless] c + c2: second radiation constant [K.cm] c + density: density of the given molecular species [molecule.cm¯³] c + drho_dP: sentivity of density to pressure [molecule.cm¯³.atm¯¹] c + drho_dT: sentivity of density to temperature [molecule.cm¯³.K¯¹] c + drho_dx: sentivity of density to the concentration [molecule.cm¯³] c + Q: total internal partition sum c + Qref: total internal partition sum at reference temperature (Tref) c + dQ_dT: sensitivity of the partition sum function to temperature c + line_profile: 1 for a Lorentz profile, 2 for a Voigt profile c + usl: options for using sub-lorentzian profiles c + sl_choice: choice of sub-lorentzian profile c + slwr: option for using CO2 SL profiles over the whole IR range c + slam: option for using CO2 SL profiles for all molecules c + slas: option for taking into account SL profile asymetry c + sens_T: option for computing sensitivities to temperature c + sens_P: option for computing sensitivities to pressure c + sens_x: option for computing sensitivities to concentrations c c Outputs: c + sigma: absorption cross-section [cm²/molecule] c + k: absorption coefficient [m¯¹] c + dsigma_dP: sensitivity of sigma to total pressure [cm²/molecule/atm] c + dsigma_dT: sensitivity of sigma to temperature [cm²/molecule/K] c + dsigma_dx: sensitivity of sigma to concentration [cm²/molecule] c + dk_dP: sensitivity of k to total pressure [m¯¹/atm] c + dk_dT: sensitivity of k to temperature [m¯¹/K] c + dk_dx: sensitivity of k to concentration [m¯¹] c include 'parameters.inc' c Inputs double precision nu double precision Tref double precision T integer index double precision P double precision Ps double precision nu0 double precision mol_mass double precision Sref double precision gamma_air double precision gamma_self double precision Elow double precision n_air double precision delta_air double precision n_self double precision c2 double precision density double precision drho_dP double precision drho_dT double precision drho_dx double precision Q double precision Qref double precision dQ_dT integer line_profile logical usl integer sl_choice logical slwr logical slam logical slas logical sens_T,sens_P,sens_x c Outputs double precision sigma,k double precision dsigma_dP,dsigma_dT,dsigma_dx double precision dk_dP,dk_dT,dk_dx c temporary double precision S,f double precision dS_dT double precision gamma_L,gamma_D double precision dgammaL_dP,dgammaL_dT,dgammaL_dx double precision dgammaD_dT integer out1,out2 double precision khi double precision nuc,dnuc_dP,x,y,k_func integer strlen character*(Nchar_mx) label label='subroutine calc_k' c Debug c write(*,*) 'from calc_k: Tref=',Tref c Debug nuc=nu0+delta_air*P dnuc_dP=delta_air c Lorentz half-width [cm¯¹] call Lorentz_linewidth(Tref,T,P,Ps,n_air,n_self, & gamma_air,gamma_self,gamma_L) c Doppler half-width [cm¯¹] call Doppler_linewidth(T,mol_mass,nu0,gamma_D) c sensitivities if ((sens_P).or.(sens_T).or.(sens_x)) then call gammaL_sensitivities(Tref,gamma_air,gamma_self, & n_air,n_self, & P,T,Ps,gamma_L, & dgammaL_dP,dgammaL_dT,dgammaL_dx) else dgammaL_dP=0.0D+0 dgammaL_dT=0.0D+0 dgammaL_dx=0.0D+0 endif if (sens_T) then call gammaD_sensitivity(mol_mass,nu0,T,gamma_D,dgammaD_dT) else dgammaD_dT=0.0D+0 endif c line intensity S [cm¯¹/(molecule.cm¯²)] call line_intensity(Tref,T,nu0,Sref,Elow,c2,Q,Qref,S) c sensitivity of line intensity to temperature dS_dT [cm¯¹/(molecule.cm¯²)/K] if (sens_T) then call S_sensitivity(Tref,T,nu0,Sref,Elow,c2,Q,dQ_dT,Qref,S, & dS_dT) else dS_dT=0.0D+0 endif call calc_k_opt(nu,nuc,dnuc_dP,P,Ps,T,index, & gamma_L,dgammaL_dP,dgammaL_dT,dgammaL_dx, & gamma_D,dgammaD_dT, & S,dS_dT, & density,drho_dP,drho_dT,drho_dx, & line_profile,usl,sl_choice,slwr,slam,slas, & sens_T,sens_P,sens_x, & sigma,k,dsigma_dP,dsigma_dT,dsigma_dx,dk_dP,dk_dT,dk_dx) c Debug c write(*,*) 'calc_k out' c Debug return end subroutine calc_k_opt(nu,nuc,dnuc_dP,P,Ps,T,index, & gamma_L,dgammaL_dP,dgammaL_dT,dgammaL_dx, & gamma_D,dgammaD_dT, & S,dS_dT, & density,drho_dP,drho_dT,drho_dx, & line_profile,usl,sl_choice,slwr,slam,slas, & sens_T,sens_P,sens_x, & sigma,k,dsigma_dP,dsigma_dT,dsigma_dx,dk_dP,dk_dT,dk_dx) implicit none include 'max.inc' c c Purpose: to compute the absorption coefficient c for a single spectral line, at a given wavenumber c c Inputs: c + nu: wavenumber [cm¯¹] c + nuc: wavenumber at the center of the line [cm¯¹] (nu0+delta_air*P) c + dnuc_dP: sensitivity of nuc to total pressure [cm¯¹/atm] c + P: pressure [atm] c + Ps: partial pressure of the given molecular species [atm] c + T: temperature (K) c + index of the molecule c + gamma_L: Lorentz half width c + dgammaL_dP: sensitivity of Lorentz line width to total pressure c + dgammaL_dT: sensitivity of Lorentz line width to temperature c + dgammaL_dx: sensitivity of Lorentz line width to concentration c + gamma_D: Doppler half width c + dgammaD_dT: sensitivity of Doppler line width to temperature c + S: line intensity [cm¯¹/(molecule.cm¯²)] c + dS_dT: sensitivity of line intensity S to temperature [cm¯¹/(molecule.cm¯²)/K] c + density: density of the molecular species c + drho_dP: sentivity of density to pressure [molecule.cm¯³.atm¯¹] c + drho_dT: sentivity of density to temperature [molecule.cm¯³.K¯¹] c + drho_dx: sentivity of density to the concentration [molecule.cm¯³] c + line_profile: 1 for a Lorentz profile, 2 for a Voigt profile c + usl: options for using sub-lorentzian profiles c + sl_choice: choice of sub-lorentzian profile c + slwr: option for using CO2 SL profiles over the whole IR range c + slam: option for using CO2 SL profiles for all molecules c + slas: option for taking into account SL profile asymetry c + sens_T: option for computing sensitivities to temperature c + sens_P: option for computing sensitivities to pressure c + sens_x: option for computing sensitivities to concentrations c c Outputs: c + sigma: absorption cross-section [cm²/molecule] that correspond to allowed transitions only c + k: absorption coefficient [m¯¹] allowed transitions c + dsigma_dP: sensitivity of sigma to total pressure [cm²/molecule/atm] c + dsigma_dT: sensitivity of sigma to temperature [cm²/molecule/K] c + dsigma_dx: sensitivity of sigma to concentration [cm²/molecule] c + dk_dP: sensitivity of k to total pressure [m¯¹/atm] c + dk_dT: sensitivity of k to temperature [m¯¹/K] c + dk_dx: sensitivity of k to concentration [m¯¹] c include 'parameters.inc' c Inputs double precision nu double precision nuc,dnuc_dP double precision P,Ps,T integer index double precision gamma_L,dgammaL_dP,dgammaL_dT,dgammaL_dx double precision gamma_D,dgammaD_dT double precision S,dS_dT double precision density double precision drho_dP double precision drho_dT double precision drho_dx integer line_profile logical usl integer sl_choice logical slwr logical slam logical slas logical sens_T,sens_P,sens_x c Outputs double precision sigma double precision k double precision dsigma_dP,dsigma_dT,dsigma_dx double precision dk_dP,dk_dT,dk_dx c temporary double precision x,y,k_func double precision f,khi double precision df_dP,df_dT,df_dx double precision T2,P2,Ps2 double precision gammaL2,gammaD2,x2,y2,k_func2,khi2,f2,nuc2 double precision dP,dT,dx integer out1,out2 integer strlen character*(Nchar_mx) label label='subroutine calc_k_opt' c c ----------------------------------------------------------------- c Select Voigt or Lorentz profile: c line_profile: 1 for Lorentz, 2 for Voigt c if (line_profile.eq.1) then c + Lorentz call Lorentz(nu,nuc,P,gamma_L,0.0D+0,f) ! [f]= [1/cm¯¹] else if (line_profile.eq.2) then c + Voigt x=dsqrt(dlog(2.0D+0))*(nu-nuc)/gamma_D y=dsqrt(dlog(2.0D+0))*gamma_L/gamma_D call Voigt(x,y,k_func) call khi_profile(nu,nuc,T,index, & usl,sl_choice,slwr,slam,slas,khi) f=k_func*dsqrt(dlog(2.0D+0)/pi)/gamma_D*khi ! [f]= [1/cm¯¹] else call error(label) write(*,*) 'line_profile=',line_profile write(*,*) 'Allowed values: 1, 2' stop endif c c ----------------------------------------------------------------- c call test_nan(f,out1) if (out1.eq.1) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'f=',f write(*,*) 'nuc=',nuc write(*,*) 'nu=',nu write(*,*) 'P=',P write(*,*) 'gamma_L=',gamma_L if (line_profile.eq.2) then write(*,*) 'k_func=',k_func write(*,*) 'gamma_D=',gamma_D write(*,*) 'khi=',khi endif stop endif c Debug c if ((nu.gt.2219.0D+0).and.(nu.lt.2220.0D+0)) then c write(*,*) nu,nuc,index,f c endif c Debug sigma=S*f ! [cm²/molecule] c -------------------------------------- c sensitivities of sigma c + to pressure if (sens_P) then dP=P/1.0D+3 ! atm P2=P+dP gammaL2=gamma_L+dgammaL_dP*dP if (line_profile.eq.1) then c + Lorentz call Lorentz(nu,nuc,P,gammaL2,0.0D+0,f2) ! [f]= [1/cm¯¹] else if (line_profile.eq.2) then c + Voigt gammaD2=gamma_D nuc2=nuc+dnuc_dP*dP x2=dsqrt(dlog(2.0D+0))*(nu-nuc2)/gammaD2 y2=dsqrt(dlog(2.0D+0))*gammaL2/gammaD2 call Voigt(x2,y2,k_func2) call khi_profile(nu,nuc2,T,index, & usl,sl_choice,slwr,slam,slas,khi2) f2=k_func2*dsqrt(dlog(2.0D+0)/pi)/gammaD2*khi2 endif ! line_profile df_dP=(f2-f)/dP dsigma_dP=S*df_dP c Debug call test_nan(dsigma_dP,out1) call test_inf(dsigma_dP,out2) if ((out1.eq.1).or.(out2.eq.1)) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'dsigma_dP=',dsigma_dP write(*,*) 'S=',S write(*,*) 'df_dP=',df_dP stop endif c Debug else dsigma_dP=0.0D+0 endif c + to temperature if (sens_T) then dT=1.0D+0 ! K T2=T+dT gammaL2=gamma_L+dgammaL_dT*dT if (line_profile.eq.1) then c + Lorentz call Lorentz(nu,nuc,P,gammaL2,0.0D+0,f2) ! [f]= [1/cm¯¹] else if (line_profile.eq.2) then c + Voigt gammaD2=gamma_D+dgammaD_dT*dT x2=dsqrt(dlog(2.0D+0))*(nu-nuc)/gammaD2 y2=dsqrt(dlog(2.0D+0))*gammaL2/gammaD2 call Voigt(x2,y2,k_func2) call khi_profile(nu,nuc,T2,index, & usl,sl_choice,slwr,slam,slas,khi2) f2=k_func2*dsqrt(dlog(2.0D+0)/pi)/gammaD2*khi2 endif ! line_profile df_dT=(f2-f)/dT dsigma_dT=dS_dT*f+S*df_dT c Debug call test_nan(dsigma_dT,out1) call test_inf(dsigma_dT,out2) if ((out1.eq.1).or.(out2.eq.1)) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'dsigma_dT=',dsigma_dT write(*,*) 'S=',S write(*,*) 'df_dT=',df_dT write(*,*) 'dS_dT=',dS_dT write(*,*) 'f=',f stop endif c Debug else dsigma_dT=0.0D+0 endif c + to concentration if (sens_x) then x=Ps/P dx=x/1.0D+3 gammaL2=gamma_L+dgammaL_dx*dx if (line_profile.eq.1) then c + Lorentz call Lorentz(nu,nuc,P,gammaL2,0.0D+0,f2) ! [f]= [1/cm¯¹] else if (line_profile.eq.2) then c + Voigt gammaD2=gamma_D x2=dsqrt(dlog(2.0D+0))*(nu-nuc)/gammaD2 y2=dsqrt(dlog(2.0D+0))*gammaL2/gammaD2 call Voigt(x2,y2,k_func2) call khi_profile(nu,nuc,T,index, & usl,sl_choice,slwr,slam,slas,khi2) f2=k_func2*dsqrt(dlog(2.0D+0)/pi)/gammaD2*khi2 endif ! line_profile df_dx=(f2-f)/dx dsigma_dx=S*df_dx c Debug call test_nan(dsigma_dx,out1) call test_inf(dsigma_dx,out2) if ((out1.eq.1).or.(out2.eq.1)) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'dsigma_dx=',dsigma_dx write(*,*) 'S=',S write(*,*) 'df_dx=',df_dx write(*,*) 'f=',f write(*,*) 'f2=',f2 write(*,*) 'dx=',dx write(*,*) 'P=',P write(*,*) 'Ps=',Ps write(*,*) 'index=',index stop endif c Debug else dsigma_dx=0.0D+0 endif c -------------------------------------- k=sigma*density ! [cm¯¹] k=k*1.0D+2 ! [m¯¹] c Debug call test_nan(k,out1) call test_inf(k,out2) if ((out1.eq.1).or.(out2.eq.1)) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'k=',k write(*,*) 'sigma=',sigma write(*,*) 'gamma_L=',gamma_L write(*,*) 'S=',S write(*,*) 'f=',f write(*,*) 'density=',density stop endif c Debug c -------------------------------------- c sensitivities of k c + to pressure if (sens_P) then dk_dP=1.0D+2*(dsigma_dP*density+sigma*drho_dP) ! [m¯¹/atm] c Debug call test_nan(dk_dP,out1) call test_inf(dk_dP,out2) if ((out1.eq.1).or.(out2.eq.1)) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'dk_dP=',dk_dP write(*,*) 'sigma=',sigma write(*,*) 'density=',density write(*,*) 'dsigma_dP=',dsigma_dP write(*,*) 'drho_dP=',drho_dP stop endif c Debug endif c + to temperature if (sens_T) then dk_dT=1.0D+2*(dsigma_dT*density+sigma*drho_dT) ! [m¯¹/K] c Debug call test_nan(dk_dT,out1) call test_inf(dk_dT,out2) if ((out1.eq.1).or.(out2.eq.1)) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'dk_dT=',dk_dT write(*,*) 'sigma=',sigma write(*,*) 'density=',density write(*,*) 'dsigma_dT=',dsigma_dT write(*,*) 'drho_dT=',drho_dT stop endif c Debug c Debug c if (dk_dT.eq.0.0D+0) then c write(*,*) 'dk_dT=',dK_dT c write(*,*) 'dsigma_dT*density=',dsigma_dT*density c write(*,*) 'sigma*drho_dT=',sigma*drho_dT c stop c endif c Debug endif c + to concentration if (sens_x) then dk_dx=1.0D+2*(dsigma_dx*density+sigma*drho_dx) ! [m¯¹] c Debug call test_nan(dk_dx,out1) call test_inf(dk_dx,out2) if ((out1.eq.1).or.(out2.eq.1)) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'dk_dx=',dk_dx write(*,*) 'sigma=',sigma write(*,*) 'density=',density write(*,*) 'dsigma_dx=',dsigma_dx write(*,*) 'drho_dx=',drho_dx stop endif c Debug endif c -------------------------------------- return end subroutine gammaL_sensitivities(Tref,gamma_air,gamma_self, & n_air,n_self, & P,T,Ps,gammaL, & dgammaL_dP,dgammaL_dT,dgammaL_dx) implicit none include 'max.inc' c c Purpose: to evaluate the sensitivities of gamma_L c to total pressure, temperature and concentration c c Input c + Tref: reference temperature c + gamma_air: air-broadened half-width [cm¯¹.atm¯¹] c + gamma_self: self-broadened half-width [cm¯¹.atm¯¹] c + n_air: temperature-dependance exponent for gamma_air [unitless] c + n_self: temperature-dependance exponent for gamma_self [unitless] c + P: total pressure c + T: temperature c + Ps: partial pressure of the species c + gammaL: Lorentz line width at P,T,Ps c c Output c + dgammaL_dP: sensitivity of Lorentz line width to total pressure c + dgammaL_dT: sensitivity of Lorentz line width to temperature c + dgammaL_dx: sensitivity of Lorentz line width to concentration c c I/O double precision Tref,gamma_air,gamma_self,n_air,n_self double precision P,T,Ps,gammaL double precision dgammaL_dP double precision dgammaL_dT double precision dgammaL_dx c temp double precision P1,P2,dP,deltaP double precision T1,T2,dT,deltaT double precision x,x1,x2,dx,deltax,Ps2 double precision gammaL1,gammaL2,dgammaL double precision diff,dmin integer niter,Niter_mx integer strlen character*(Nchar_mx) label label='subroutine Q_sensitivity' dmin=1.0D-3 Niter_mx=10 P1=P deltaP=P/1.0D+3 ! atm T1=T deltaT=1.0D+0 ! K if (P.ne.0.0D+0) then x=Ps/P else write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'P=',P stop endif x1=x deltax=x/1.0D+3 gammaL1=gammaL c sensitivity to pressure diff=0.0D+0 niter=0 do while (diff.lt.dmin) niter=niter+1 P2=P1+deltaP*niter dP=P2-P1 call Lorentz_linewidth(Tref,T,P2,Ps,n_air,n_self, & gamma_air,gamma_self,gammaL2) dgammaL=gammaL2-gammaL1 diff=dabs(dgammaL/gammaL1) if (niter.gt.Niter_mx) then c write(*,*) 'Error from ',label(1:strlen(label)),' :' c write(*,*) 'Maximum number of iterations reached' c write(*,*) 'diff=',diff c stop goto 121 endif enddo 121 continue dgammaL_dP=dgammaL/dP c sensitivity to temperature diff=0.0D+0 niter=0 do while (diff.lt.dmin) niter=niter+1 T2=T1+deltaT*niter dT=T2-T1 call Lorentz_linewidth(Tref,T2,P,Ps,n_air,n_self, & gamma_air,gamma_self,gammaL2) dgammaL=gammaL2-gammaL1 diff=dabs(dgammaL/gammaL1) if (niter.gt.Niter_mx) then c write(*,*) 'Error from ',label(1:strlen(label)),' :' c write(*,*) 'Maximum number of iterations reached' c write(*,*) 'diff=',diff c stop goto 122 endif enddo 122 continue dgammaL_dT=dgammaL/dT c Debug c write(*,*) 'ir dgammaL_dT=',dgammaL_dT c Debug c sensitivity to concentration diff=0.0D+0 niter=0 do while (diff.lt.dmin) niter=niter+1 x2=x1+deltax*niter dx=x2-x1 Ps2=x2*P call Lorentz_linewidth(Tref,T,P,Ps2,n_air,n_self, & gamma_air,gamma_self,gammaL2) dgammaL=gammaL2-gammaL1 diff=dabs(dgammaL/gammaL1) if (niter.gt.Niter_mx) then c write(*,*) 'Error from ',label(1:strlen(label)),' :' c write(*,*) 'Maximum number of iterations reached' c write(*,*) 'diff=',diff c stop goto 123 endif enddo 123 continue dgammaL_dx=dgammaL/dx return end subroutine gammaD_sensitivity(mol_mass,nu0,T,gammaD,dgammaD_dT) implicit none include 'max.inc' c c Purpose: to evaluate the sensitivity of gamma_D c to temperature c c Input c + mol_mass: molar mass of the molecular species [g/mol] c + nu0: vacuum wavenumber [cm¯¹] c + T: temperature c + gammaD: Doppler line width at T c c Output c + dgammaD_dT: sensitivity of Doppler line width to temperature c c I/O double precision mol_mass,nu0 double precision T,gammaD double precision dgammaD_dT c temp double precision T1,T2,dT,deltaT double precision gammaD1,gammaD2,dgammaD double precision diff,dmin integer niter,Niter_mx integer strlen character*(Nchar_mx) label label='subroutine Q_sensitivity' dmin=1.0D-4 Niter_mx=10 T1=T deltaT=1.0D+0 ! K gammaD1=gammaD diff=0.0D+0 niter=0 do while (diff.lt.dmin) niter=niter+1 T2=T1+deltaT*niter dT=T2-T1 call Doppler_linewidth(T2,mol_mass,nu0,gammaD2) dgammaD=gammaD2-gammaD1 diff=dabs(dgammaD/gammaD1) if (niter.gt.Niter_mx) then c write(*,*) 'Error from ',label(1:strlen(label)),' :' c write(*,*) 'Maximum number of iterations reached' c write(*,*) 'diff=',diff c stop goto 123 endif enddo 123 continue dgammaD_dT=dgammaD/dT c Debug c write(*,*) 'ir dgammaD_dT=',dgammaD_dT c Debug return end subroutine Q_sensitivity(index,isotope,iso_g,T,Q,dQ_dT) implicit none include 'max.inc' c c Purpose: to evaluate the sensitivity of partition function Q c to temperature c c Input c + index: number of the HITRAN molecule c + isotope: number of the HITRAN isotope c + iso_g: state independent degeneracy factor c + T: temperature c + Q: partition function at T c c Output c + dQ_dT: sensitivity of Q to temperature c c I/O integer index integer isotope integer iso_g double precision T,Q double precision dQ_dT c temp double precision T1,T2,dT,deltaT double precision Q1,Q2,dQ double precision diff,dmin integer niter,Niter_mx integer strlen character*(Nchar_mx) label label='subroutine Q_sensitivity' dmin=1.0D-4 Niter_mx=10 T1=T deltaT=1.0D+0 ! K Q1=Q diff=0.0D+0 niter=0 do while (diff.lt.dmin) niter=niter+1 T2=T1+deltaT*niter dT=T2-T1 call TIPS(index,T2,isotope,iso_g,Q2) dQ=Q2-Q1 diff=dabs(dQ/Q1) if (niter.gt.Niter_mx) then c write(*,*) 'Error from ',label(1:strlen(label)),' :' c write(*,*) 'Maximum number of iterations reached' c write(*,*) 'diff=',diff c stop goto 123 endif enddo 123 continue dQ_dT=dQ/dT c Debug c write(*,*) 'ir dQ_dT=',dQ_dT c Debug return end subroutine S_sensitivity(Tref,T,nu0,Sref,Elow,c2,Q,dQ_dT,Qref,S, & dS_dT) implicit none include 'max.inc' c c Purpose: to evaluate the sensitivities of line intensity S c to temperature. c c Input c + Tref: reference temperature [K] c + T: temperature [K] c + nu0: vacuum wavenumber [cm¯¹] c + Sref: line intensity [cm¯¹/molecule.cm¯²] at standard Tref c + Elow: lower-state energy [cm¯¹] c + c2: second radiation constant [K.cm] c + Q: total internal partition sum at T c + dQ_dT: sensitivity of Q to T c + Qref: total internal partition sum at reference temperature (Tref) c + S: line intensity at temperature T [cm¯¹/(molecule.cm¯²)] c c Output c + dS_dT: sensitivity of S to temperature [cm¯¹/(molecule.cm¯²)/K] c c I/O double precision Tref,T,nu0,Sref,Elow,c2,Q,dQ_dT,Qref,S double precision dS_dT c temp double precision T1,T2,dT,deltaT,Q2 double precision S1,S2,dS double precision diff,dmin integer niter,Niter_mx integer strlen character*(Nchar_mx) label label='subroutine S_sensitivities' dmin=1.0D-4 Niter_mx=10 T1=T deltaT=1.0D+0 ! K S1=S diff=0.0D+0 niter=0 do while (diff.lt.dmin) niter=niter+1 T2=T1+deltaT*niter dT=T2-T1 Q2=Q+dQ_dT*dT call line_intensity(Tref,T2,nu0,Sref,Elow,c2,Q2,Qref,S2) dS=S2-S1 diff=dabs(dS/S1) if (niter.gt.Niter_mx) then c write(*,*) 'Error from ',label(1:strlen(label)),' :' c write(*,*) 'Maximum number of iterations reached' c write(*,*) 'diff=',diff c stop goto 123 endif enddo 123 continue dS_dT=dS/dT c Debug c write(*,*) 'ir dS_dT=',dS_dT c Debug return end