subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall,ind) !================================================================== ! ! Purpose ! ------- ! Calculates the H2-He CIA opacity, using a lookup table from ! HITRAN (2011) ! ! Authors ! ------- ! R. Wordsworth (2011) ! !================================================================== use callkeys_mod, only: H2orthopara_mixture, strictboundcia use datafile_mod, only: datadir use mod_phys_lmdz_para, only : is_master 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(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 temp_arr(nT) double precision, dimension(:), allocatable :: wn_arr double precision, dimension(:,:), allocatable :: abs_arr integer k,iT logical firstcall save ns, wn_arr, temp_arr, abs_arr !read by master 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, ztemp integer nres integer ind ztemp=temp if (H2orthopara_mixture .eq. "hot") then ! .and. (temp .gt. 3000.0)) then ! print*,"We're in the Hot Jupiter case " if (ztemp .gt. 9900) then if (strictboundcia) then if (is_master) then print*,'Your temperatures are too high for this H2-He CIA dataset (Hot Jupiter case).' print*,'Please run mixed H2-He atmospheres below T = 9900.0 K.' endif !is_master stop else !if (is_master) then ! print*,'Your temperatures are too high for this H2-He CIA dataset (Hot Jupiter case)' ! print*,'you have chosen strictboundcia = ', strictboundcia ! print*,'*********************************************************' ! print*,' we allow model to continue but with temp = 9900.0 K ' ! print*,' ... for H2-He CIA dataset ... ' ! print*,' ... we assume we know what you are doing ... ' ! print*,'*********************************************************' !endif !is_master ztemp = 9900.0 endif ! of stricbound cia endif ! temp .gt. 9900 else !not in Hot Jupiter case if(ztemp.gt.400)then if (strictboundcia) then if (is_master) then print*,'Your temperatures are too high for this H2-He CIA dataset.' print*,'Please run mixed H2-He atmospheres below T = 400 K.' endif ! is_master stop else !if (is_master) then ! print*,'Your temperatures are too high for this H2-He CIA dataset' ! print*,'you have chosen strictboundcia = ', strictboundcia ! print*,'*********************************************************' ! print*,' we allow model to continue but with temp = 400 ' ! print*,' ... for H2-He CIA dataset ... ' ! print*,' ... we assume we know what you are doing ... ' ! print*,'*********************************************************' !endif ! is_master ztemp = 400 endif !of strictboundcia elseif(ztemp.lt.40)then if (strictboundcia) then if (is_master) then print*,'Your temperatures are too low for this H2-He CIA dataset.' print*,'Please run mixed H2-He atmospheres above T = 40 K.' endif ! is_master stop else !if (is_master) then ! print*,'Your temperatures are too low for this H2-He CIA dataset' ! print*,'you have chosen strictboundcia = ', strictboundcia ! print*,'*********************************************************' ! print*,' we allow model to continue but with temp = 40 ' ! print*,' ... for H2-He CIA dataset ... ' ! print*,' ... we assume we know what you are doing ... ' ! print*,'*********************************************************' !endif ! is_master ztemp = 40 endif !of strictboundcia endif ! of (temp .gt. 400) endif ! of (H2orthopara_mixture .eq. "hot") amagatH2 = (273.15/temp)*(presH2/101325.0) amagatHe = (273.15/temp)*(presHe/101325.0) if(firstcall)then ! called by sugas_corrk only if (is_master) print*,'----------------------------------------------------' if (is_master) print*,'Initialising H2-He continuum from HITRAN database...' ! 1.1 Open the ASCII files if (H2orthopara_mixture.eq."normal") then nS = 2428 dt_file=TRIM(datadir)//'/continuum_data/H2-He_norm_2011.cia' else if (H2orthopara_mixture.eq."equilibrium") then nS = 1501 dt_file=TRIM(datadir)//'/continuum_data/H2-He_eq_0-15000cm-1_40-400K.cia' else if (H2orthopara_mixture .eq. "hot") then !we use a dataset than can go as high as 9900 K nS = 19981 dt_file=TRIM(datadir)//'/continuum_data/H2-He_2011.cia' endif if (.not. allocated(wn_arr)) allocate(wn_arr(nS)) if (.not. allocated(abs_arr)) allocate(abs_arr(nS,nT)) !$OMP MASTER open(33,file=dt_file,form='formatted',status='old',iostat=ios) if (ios.ne.0) then ! file not found if (is_master) then 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 is there.' endif call abort else do iT=1,nT read(33,fmat1) bleh,blah,blah,nres,Ttemp if(nS.ne.nres)then if (is_master) 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' endif 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) !$OMP END MASTER !$OMP BARRIER if (is_master) then print*,'interpolateH2He: At wavenumber ',wn,' cm^-1' print*,' temperature ',temp,' K' print*,' H2 partial pressure ',presH2,' Pa' print*,' and He partial pressure ',presHe,' Pa' endif endif !! if (firstcall) call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,ztemp,abcoef,ind) !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' ! 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. return end subroutine interpolateH2He