subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall) !================================================================== ! ! Purpose ! ------- ! Calculates the H2-He CIA opacity, using a lookup table from ! HITRAN (2011) ! ! Authors ! ------- ! R. Wordsworth (2011) ! !================================================================== use datafile_mod, only: datadir implicit none ! input double precision wn ! wavenumber (cm^-1) double precision temp ! temperature (Kelvin) double precision presH2 ! H2 partial pressure (Pascals) double precision presHe ! He partial pressure (Pascals) ! output double precision abcoef ! absorption coefficient (m^-1) integer nS,nT parameter(nS=2428) parameter(nT=10) double precision, parameter :: losch = 2.6867774e19 ! Loschmit's number (molecule cm^-3 at STP) ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2 ! see Richard et al. 2011, JQSRT for details double precision amagatH2 double precision amagatHe double precision wn_arr(nS) double precision temp_arr(nT) double precision abs_arr(nS,nT) integer k,iT logical firstcall save wn_arr, temp_arr, abs_arr character*100 dt_file integer strlen,ios character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" character*20 bleh double precision blah, Ttemp integer nres if(temp.gt.400)then print*,'Your temperatures are too high for this H2-He CIA dataset.' print*,'Please run mixed H2-He atmospheres below T = 400 K.' stop endif amagatH2 = (273.15/temp)*(presH2/101325.0) amagatHe = (273.15/temp)*(presHe/101325.0) if(firstcall)then ! called by sugas_corrk only print*,'----------------------------------------------------' print*,'Initialising H2-He continuum from HITRAN database...' ! 1.1 Open the ASCII files dt_file=TRIM(datadir)//'/continuum_data/H2-He_norm_2011.cia' open(33,file=dt_file,form='formatted',status='old',iostat=ios) if (ios.ne.0) then ! file not found write(*,*) 'Error from interpolateH2He' write(*,*) 'Data file ',trim(dt_file),' not found.' write(*,*) 'Check that your path to datagcm:',trim(datadir) write(*,*) 'is correct. You can change it in callphys.def with:' write(*,*) 'datadir = /absolute/path/to/datagcm' write(*,*) 'Also check that the continuum data continuum_data/H2-He_norm_2011.cia is there.' call abort else do iT=1,nT read(33,fmat1) bleh,blah,blah,nres,Ttemp if(nS.ne.nres)then print*,'Resolution given in file: ',trim(dt_file) print*,'is ',nres,', which does not match nS.' print*,'Please adjust nS value in interpolateH2He.F90' stop endif temp_arr(iT)=Ttemp do k=1,nS read(33,*) wn_arr(k),abs_arr(k,it) end do end do endif close(33) print*,'interpolateH2He: At wavenumber ',wn,' cm^-1' print*,' temperature ',temp,' K' print*,' H2 partial pressure ',presH2,' Pa' print*,' and He partial pressure ',presHe,' Pa' call bilinearH2He(wn_arr,temp_arr,abs_arr,wn,temp,abcoef) print*,'the absorption is ',abcoef,' cm^5 molecule^-2' print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' abcoef=abcoef*losch**2*100.0*amagatH2*amagatHe ! convert to m^-1 print*,'We have ',amagatH2,' amagats of H2' print*,'and ',amagatHe,' amagats of He' print*,'So the absorption is ',abcoef,' m^-1' else call bilinearH2He(wn_arr,temp_arr,abs_arr,wn,temp,abcoef) abcoef=abcoef*losch**2*100.0*amagatH2*amagatHe ! convert to m^-1 ! unlike for Rayleigh scattering, we do not currently weight by the BB function ! however our bands are normally thin, so this is no big deal. endif return end subroutine interpolateH2He !------------------------------------------------------------------------- subroutine bilinearH2He(x_arr,y_arr,f2d_arr,x_in,y_in,f) ! Necessary for interpolation of continuum data implicit none integer nX,nY,i,j,a,b parameter(nX=2428) parameter(nY=10) real*8 x_in,y_in,x,y,x1,x2,y1,y2 real*8 f,f11,f12,f21,f22,fA,fB real*8 x_arr(nX) real*8 y_arr(nY) real*8 f2d_arr(nX,nY) integer strlen character*100 label label='subroutine bilinear' x=x_in y=y_in ! 1st check we're within the wavenumber range if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then f=0.0D+0 return else ! in the x (wavenumber) direction 1st i=1 10 if (i.lt.(nX+1)) then if (x_arr(i).gt.x) then x1=x_arr(i-1) x2=x_arr(i) a=i-1 i=9999 endif i=i+1 goto 10 endif endif if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then write(*,*) 'Warning from bilinearH2He:' write(*,*) 'Outside continuum temperature range!' if(y.lt.y_arr(1))then y=y_arr(1)+0.01 endif if(y.gt.y_arr(nY))then y=y_arr(nY)-0.01 endif else ! in the y (temperature) direction 2nd j=1 20 if (j.lt.(nY+1)) then if (y_arr(j).gt.y) then y1=y_arr(j-1) y2=y_arr(j) b=j-1 j=9999 endif j=j+1 goto 20 endif endif f11=f2d_arr(a,b) f21=f2d_arr(a+1,b) f12=f2d_arr(a,b+1) f22=f2d_arr(a+1,b+1) call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) return end subroutine bilinearH2He