Index: trunk/LMDZ.GENERIC/changelog.txt =================================================================== --- trunk/LMDZ.GENERIC/changelog.txt (revision 4054) +++ trunk/LMDZ.GENERIC/changelog.txt (revision 4055) @@ -2168,2 +2168,6 @@ the black body function. +== 06/02/2026 == MT +Make the 'new generic continuum database' the standard and only way to handle continuum absorption. +Remove all the interpolateXY.F90 routines. Remove the bilinear and bilinearbig routines. RIP. +More details here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Continuum_Database Index: trunk/LMDZ.GENERIC/deftank/callphys.GJ581d =================================================================== --- trunk/LMDZ.GENERIC/deftank/callphys.GJ581d (revision 4054) +++ trunk/LMDZ.GENERIC/deftank/callphys.GJ581d (revision 4055) @@ -32,5 +32,4 @@ # Include continuum absorption in radiative transfer (note CO2 is treated separately) ? continuum = .true. -generic_continuum_database = .true. # folder in which correlated-k data is stored ? corrkdir = CO2_H2Ovar Index: trunk/LMDZ.GENERIC/deftank/callphys.earlymars =================================================================== --- trunk/LMDZ.GENERIC/deftank/callphys.earlymars (revision 4054) +++ trunk/LMDZ.GENERIC/deftank/callphys.earlymars (revision 4055) @@ -32,5 +32,4 @@ # Include continuum absorption in radiative transfer (note CO2 is treated separately) ? continuum = .true. -generic_continuum_database = .true. # folder in which correlated-k data is stored ? corrkdir = earlymars Index: trunk/LMDZ.GENERIC/deftank/callphys.earthslab =================================================================== --- trunk/LMDZ.GENERIC/deftank/callphys.earthslab (revision 4054) +++ trunk/LMDZ.GENERIC/deftank/callphys.earthslab (revision 4055) @@ -51,5 +51,4 @@ # call continuum in radiative transfer ? continuum = .true. -generic_continuum_database = .true. # Include Rayleigh scattering in the visible ? rayleigh = .true. Index: trunk/LMDZ.GENERIC/deftank/callphys.kcm1d =================================================================== --- trunk/LMDZ.GENERIC/deftank/callphys.kcm1d (revision 4054) +++ trunk/LMDZ.GENERIC/deftank/callphys.kcm1d (revision 4055) @@ -16,5 +16,4 @@ # Include continuum absorption in radiative transfer (note CO2 is treated separately) ? continuum = .true. -generic_continuum_database = .true. # folder in which correlated-k data is stored ? corrkdir = N2_CO2step_H2Ovar/1E3ppmCO2 Index: trunk/LMDZ.GENERIC/libf/phygeneric/bilinear.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/bilinear.F90 (revision 4054) +++ (revision ) @@ -1,18 +1,0 @@ -!------------------------------------------------------------------------- -subroutine bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) -! Used for interpolation of continuum data - - implicit none - - real*8 x,y,x1,x2,y1,y2 - real*8 f,f11,f12,f21,f22,fA,fB - - ! 1st in x-direction - fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1) - fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1) - - ! then in y-direction - f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1) - - return -end subroutine bilinear Index: trunk/LMDZ.GENERIC/libf/phygeneric/bilinearbig.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/bilinearbig.F90 (revision 4054) +++ (revision ) @@ -1,92 +1,0 @@ - subroutine bilinearbig(nX,nY,x_arr,y_arr,f2d_arr,x_in,y_in,f,ind) - -! Necessary for interpolation of continuum data -! optimized by A. Spiga 01/2013 - - implicit none - - integer nX,nY,i,j,ind,b - - real*8 x_in,y_in,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) - real*8,save :: x,y -!$OMP THREADPRIVATE(x,y) - - integer strlen - character*100 label - label='subroutine bilinear' - - - x=x_in - y=y_in - - !! AS: important to optimize here because the array is quite large - !! ... and actually calculations only need to be done once - !! IF ind=-9999 we have not calculated yet - if ( ind == -9999) then - !1st check we're within the wavenumber range - if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then - ind=-1 - else - i=1 - x2=x_arr(i) - do while ( x2 .le. x ) - x1=x2 - i=i+1 - x2=x_arr(i) - ind=i-1 - end do - endif - endif - - !! Either we already saw we are out of wavenumber range - !! ... and we just have to set f=0 and exit - if ( ind == -1) then - f=0.0D+0 - return - !! Or we already determined ind -- so we just proceed - else - x1=x_arr(ind) - x2=x_arr(ind+1) - endif - -! ... and for y within the temperature range - if ((y.le.y_arr(1)).or.(y.ge.y_arr(nY))) then - !print*,y_arr(1),y_arr(nY) - !write(*,*) 'Warning from bilinearbig routine:' - !write(*,*) 'Outside continuum temperature range!' - if(y.le.y_arr(1))then - y=y_arr(1)+0.01 - b=1 - y1=y_arr(b) - y2=y_arr(b+1) - endif - if(y.ge.y_arr(nY))then - y=y_arr(nY)-0.01 - b=nY-1 - y1=y_arr(b) - y2=y_arr(b+1) - endif - else - j=1 - y2=y_arr(j) - do while ( y2 .lt. y ) - y1=y2 - j=j+1 - y2=y_arr(j) - b=j-1 - end do - endif - - f11=f2d_arr(ind,b) - f21=f2d_arr(ind+1,b) - f12=f2d_arr(ind,b+1) - f22=f2d_arr(ind+1,b+1) - - call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) - - return - end subroutine bilinearbig Index: trunk/LMDZ.GENERIC/libf/phygeneric/callkeys_mod.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/callkeys_mod.F90 (revision 4054) +++ trunk/LMDZ.GENERIC/libf/phygeneric/callkeys_mod.F90 (revision 4055) @@ -13,6 +13,6 @@ logical,save :: season,diurnal,tlocked,rings_shadow,lwrite !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite) - logical,save :: callgasvis,continuum,generic_continuum_database,graybody -!$OMP THREADPRIVATE(callgasvis,continuum,generic_continuum_database,graybody) + logical,save :: callgasvis,continuum,graybody +!$OMP THREADPRIVATE(callgasvis,continuum,graybody) logical,save :: strictboundcorrk !$OMP THREADPRIVATE(strictboundcorrk) Index: trunk/LMDZ.GENERIC/libf/phygeneric/inifis_mod.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/inifis_mod.F90 (revision 4054) +++ trunk/LMDZ.GENERIC/libf/phygeneric/inifis_mod.F90 (revision 4055) @@ -381,12 +381,4 @@ call getin_p("continuum",continuum) if (is_master) write(*,*) trim(rname)//": continuum = ",continuum - - if (is_master) write(*,*) trim(rname)//& - ": use generic continuum opacity database ?"//& - " (matters only if callrad=T)" - generic_continuum_database=.true. ! default value - call getin_p("generic_continuum_database",generic_continuum_database) - if (is_master) write(*,*) trim(rname)//": generic_continuum_database = ", & - generic_continuum_database if (is_master) write(*,*) trim(rname)//": version for H2H2 CIA file ?" Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateCH4CH4.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateCH4CH4.F90 (revision 4054) +++ (revision ) @@ -1,170 +1,0 @@ -subroutine interpolateCH4CH4(wn,temp,pres,abcoef,firstcall,ind) - - !================================================================== - ! - ! Purpose - ! ------- - ! Calculates the CH4-CH4 CIA opacity, using a lookup table from - ! HITRAN (2011) - ! - ! Authors - ! ------- - ! R. Wordsworth (2011) - ! - !================================================================== - - use callkeys_mod, only: 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 pres ! CH4 partial pressure (Pascals) - - ! output - double precision abcoef ! absorption coefficient (m^-1) - - integer nS,nT - parameter(nS=1018) - 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 amagat - 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 !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(ztemp.gt.400)then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too high for this CH4-CH4 CIA dataset.' - print*,'Please run mixed CH4-CH4 atmospheres below T = 400 K.' - endif - call abort_physic("interpolateCH4CH4", "temperatures are too high for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too high for this CH4-CH4 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 400 ' - ! print*,' ... for CH4-CH4 CIA dataset ... ' - ! print*,' ... we assume we know what you are doing ... ' - ! print*,'*********************************************************' - !endif - ztemp = 400 - endif - elseif(ztemp.lt.40)then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too low for this CH4-CH4 CIA dataset.' - print*,'Please run mixed CH4-CH4 atmospheres above T = 40 K.' - endif - call abort_physic("interpolateCH4CH4", "temperatures are too low for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too low for this CH4-CH4 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 40 ' - ! print*,' ... for CH4-CH4 CIA dataset ... ' - ! print*,' ... we assume we know what you are doing ... ' - ! print*,'*********************************************************' - !endif - ztemp = 40 - endif - endif - - amagat = (273.15/temp)*(pres/101325.0) - - if(firstcall)then ! called by sugas_corrk only - if (is_master) print*,'----------------------------------------------------' - if (is_master) print*,'Initialising CH4-CH4 continuum from HITRAN database...' - - ! 1.1 Open the ASCII files - dt_file=TRIM(datadir)//'/continuum_data/CH4-CH4_2011.cia' - -!$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 interpolateCH4CH4' - 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/CH4-CH4_2011.cia is there.' - endif - call abort_physic("interpolateCH4CH4", "Unable to read file", 1) - 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 interpolateCH4CH4.F90' - endif - call abort_physic("interpolateCH4CH4", "Resolution given does not match with nS value", 1) - 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*,'interpolateCH4CH4: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' pressure ',pres,' Pa' - endif - endif - - 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*amagat**2 ! convert to m^-1 - -! print*,'We have ',amagat,' amagats of CH4' -! 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 interpolateCH4CH4 - Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateCO2CH4.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateCO2CH4.F90 (revision 4054) +++ (revision ) @@ -1,140 +1,0 @@ -subroutine interpolateCO2CH4(wn,temp,presCO2,presCH4,abcoef,firstcall,ind) - - !================================================================== - ! - ! Purpose - ! ------- - ! Calculates the CO2-CH4 CIA opacity, using a lookup table from - ! Turbet et al. 2020, Icarus, Volume 346, article id. 113762 - ! - ! Authors - ! ------- - ! M. Turbet (2023) - ! - !================================================================== - - use datafile_mod, only: datadir - implicit none - - ! input - double precision wn ! wavenumber (cm^-1) - double precision temp ! temperature (Kelvin) - double precision presCO2 ! CO2 partial pressure (Pascals) - double precision presCH4 ! CH4 partial pressure (Pascals) - - ! output - double precision abcoef ! absorption coefficient (m^-1) - - integer nS,nT - parameter(nS=240) - parameter(nT=6) - 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 amagatCO2 - double precision amagatCH4 - 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 !read by master - - character*100 dt_file - integer strlen,ios - - character(LEN=*), parameter :: fmat1 = "(A21,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" - - character*21 bleh - double precision blah, Ttemp - integer nres - integer ind - - if(temp.gt.600)then - print*,'Your temperatures are too high for this CO2-CH4 CIA dataset.' - print*,'Please run mixed CO2-CH4 atmospheres below T = 600 K.' - call abort_physic("interpolateCO2CH4", "temperatures are too high for this CIA dataset", 1) - endif - - if(temp.lt.100)then - print*,'Your temperatures are too low for this CO2-CH4 CIA dataset.' - print*,'Please run mixed CO2-CH4 atmospheres above T = 100 K.' - call abort_physic("interpolateCO2CH4", "temperatures are too low for this CIA dataset", 1) - endif - - amagatCO2 = (273.15/temp)*(presCO2/101325.0) - amagatCH4 = (273.15/temp)*(presCH4/101325.0) - - if(firstcall)then ! called by sugas_corrk only - print*,'----------------------------------------------------' - print*,'Initialising CO2-CH4 continuum from Turbet et al. 2020' - - ! 1.1 Open the ASCII files - dt_file=TRIM(datadir)//'/continuum_data/CO2-CH4_TURBET2020.cia' - - -!$OMP MASTER - open(33,file=dt_file,form='formatted',status='old',iostat=ios) - if (ios.ne.0) then ! file not found - write(*,*) 'Error from interpolateCO2CH4' - 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/CO2-CH4_TURBET2020.cia is there.' - call abort_physic("interpolateCO2CH4", "Unable to read file", 1) - 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 interpolateCO2CH4.F90' - call abort_physic("interpolateCO2CH4", "Resolution given does not match with nS value", 1) - endif - temp_arr(iT)=Ttemp - - do k=1,nS - read(33,*) wn_arr(k),abs_arr(k,it) - write(*,*) 'for k = ', k, ' we have ', wn_arr(k),abs_arr(k,it) - end do - - end do - - endif - close(33) -!$OMP END MASTER -!$OMP BARRIER - - print*,'interpolateCO2CH4: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' CO2 partial pressure ',presCO2,' Pa' - print*,' and CH4 partial pressure ',presCH4,' Pa' - - endif - - call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) - -! print*,'the absorption is ',abcoef,' cm^5 molecule^-2' -! print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' - - abcoef=abcoef*100.0*amagatCO2*amagatCH4 ! convert to m^-1 - -! print*,'We have ',amagatCO2,' amagats of CO2' -! print*,'and ',amagaCH4,' amagats of CH4' -! 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 interpolateCO2CH4 - Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateCO2H2.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateCO2H2.F90 (revision 4054) +++ (revision ) @@ -1,140 +1,0 @@ -subroutine interpolateCO2H2(wn,temp,presCO2,presH2,abcoef,firstcall,ind) - - !================================================================== - ! - ! Purpose - ! ------- - ! Calculates the CO2-H2 CIA opacity, using a lookup table from - ! Turbet et al. 2020, Icarus, Volume 346, article id. 113762 - ! - ! Authors - ! ------- - ! M. Turbet (2023) - ! - !================================================================== - - use datafile_mod, only: datadir - implicit none - - ! input - double precision wn ! wavenumber (cm^-1) - double precision temp ! temperature (Kelvin) - double precision presCO2 ! CO2 partial pressure (Pascals) - double precision presH2 ! H2 partial pressure (Pascals) - - ! output - double precision abcoef ! absorption coefficient (m^-1) - - integer nS,nT - parameter(nS=300) - parameter(nT=6) - 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 amagatCO2 - double precision amagatH2 - 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 !read by master - - character*100 dt_file - integer strlen,ios - - character(LEN=*), parameter :: fmat1 = "(A21,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" - - character*21 bleh - double precision blah, Ttemp - integer nres - integer ind - - if(temp.gt.600)then - print*,'Your temperatures are too high for this CO2-H2 CIA dataset.' - print*,'Please run mixed CO2-H2 atmospheres below T = 600 K.' - call abort_physic("interpolateCO2H2", "temperatures are too high for this CIA dataset", 1) - endif - - if(temp.lt.100)then - print*,'Your temperatures are too low for this CO2-H2 CIA dataset.' - print*,'Please run mixed CO2-H2 atmospheres above T = 100 K.' - call abort_physic("interpolateCO2H2", "temperatures are too low for this CIA dataset", 1) - endif - - amagatCO2 = (273.15/temp)*(presCO2/101325.0) - amagatH2 = (273.15/temp)*(presH2/101325.0) - - if(firstcall)then ! called by sugas_corrk only - print*,'----------------------------------------------------' - print*,'Initialising CO2-H2 continuum from Turbet et al. 2020' - - ! 1.1 Open the ASCII files - dt_file=TRIM(datadir)//'/continuum_data/CO2-H2_TURBET2020.cia' - - -!$OMP MASTER - open(33,file=dt_file,form='formatted',status='old',iostat=ios) - if (ios.ne.0) then ! file not found - write(*,*) 'Error from interpolateCO2H2' - 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/CO2-H2_TURBET2020.cia is there.' - call abort_physic("interpolateCO2H2", "Unable to read file", 1) - 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 interpolateCO2H2.F90' - call abort_physic("interpolateCO2H2", "Resolution given does not match with nS value", 1) - endif - temp_arr(iT)=Ttemp - - do k=1,nS - read(33,*) wn_arr(k),abs_arr(k,it) - write(*,*) 'for k = ', k, ' we have ', wn_arr(k),abs_arr(k,it) - end do - - end do - - endif - close(33) -!$OMP END MASTER -!$OMP BARRIER - - print*,'interpolateCO2H2: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' CO2 partial pressure ',presCO2,' Pa' - print*,' and H2 partial pressure ',presH2,' Pa' - - endif - - call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind) - -! print*,'the absorption is ',abcoef,' cm^5 molecule^-2' -! print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2' - - abcoef=abcoef*100.0*amagatCO2*amagatH2 ! convert to m^-1 - -! print*,'We have ',amagatCO2,' amagats of CO2' -! print*,'and ',amagatH2,' amagats of H2' -! 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 interpolateCO2H2 - Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateH2CH4.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateH2CH4.F90 (revision 4054) +++ (revision ) @@ -1,179 +1,0 @@ -subroutine interpolateH2CH4(wn,temp,presH2,presCH4,abcoef,firstcall,ind) - - !================================================================== - ! - ! Purpose - ! ------- - ! Calculates the H2-CH4 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 presCH4 ! CH4 partial pressure (Pascals) - double precision presH2 ! H2 partial pressure (Pascals) - - ! output - double precision abcoef ! absorption coefficient (m^-1) - - integer nS,nT - parameter(nS=1974) - 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 amagatCH4 - double precision amagatH2 - 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 !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(ztemp.gt.400)then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too high for this H2-CH4 CIA dataset.' - print*,'Please run mixed H2-CH4 atmospheres below T = 400 K.' - endif - call abort_physic("interpolateH2CH4", "temperatures are too high for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too high for this H2-CH4 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 400 ' - ! print*,' ... for H2-CH4 CIA dataset ... ' - ! print*,' ... we assume we know what you are doing ... ' - ! print*,'*********************************************************' - !endif - ztemp = 400 - endif - elseif(ztemp.lt.40)then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too low for this H2-CH4 CIA dataset.' - print*,'Please run mixed H2-CH4 atmospheres above T = 40 K.' - endif - call abort_physic("interpolateH2CH4", "temperatures are too low for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too low for this H2-CH4 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 40 ' - ! print*,' ... for H2-CH4 CIA dataset ... ' - ! print*,' ... we assume we know what you are doing ... ' - ! print*,'*********************************************************' - !endif - ztemp = 40 - endif - endif - - amagatCH4 = (273.15/temp)*(presCH4/101325.0) - amagatH2 = (273.15/temp)*(presH2/101325.0) - - if(firstcall)then ! called by sugas_corrk only - if (is_master) print*,'----------------------------------------------------' - if (is_master) print*,'Initialising H2-CH4 continuum from HITRAN database...' - - ! 1.1 Open the ASCII files - if (H2orthopara_mixture.eq."normal") then - dt_file=TRIM(datadir)//'/continuum_data/H2-CH4_norm_2011.cia' - else if ((H2orthopara_mixture.eq."hot") .or. (H2orthopara_mixture.eq."equilibrium")) then - dt_file=TRIM(datadir)//'/continuum_data/H2-CH4_eq_2011.cia' - endif - -!$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 interpolateH2CH4' - 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-CH4_norm_2011.cia is there.' - endif - call abort_physic("interpolateH2CH4", "Unable to read file", 1) - 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 interpolateH2CH4.F90' - endif - call abort_physic("interpolateH2CH4", "Resolution given does not match with nS value", 1) - 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*,'interpolateH2CH4: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' CH4 partial pressure ',presCH4,' Pa' - print*,' and H2 partial pressure ',presH2,' Pa' - endif - endif - - 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*amagatCH4*amagatH2 ! convert to m^-1 - -! print*,'We have ',amagatCH4,' amagats of CH4' -! print*,'and ',amagatH2,' amagats of H2' -! 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 interpolateH2CH4 - Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateH2H2.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateH2H2.F90 (revision 4054) +++ (revision ) @@ -1,234 +1,0 @@ - - subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,ind) - -!================================================================== -! -! Purpose -! ------- -! Calculates the H2-H2 CIA opacity, using a lookup table from -! HITRAN (2011 or later) -! -! Authors -! ------- -! R. Wordsworth (2011) -! -! + J.Vatant d'Ollone (2019) -! - Enable updated HITRAN file (Karman2019,Fletcher2018) -! (2018 one should be default for giant planets) -!================================================================== - - use callkeys_mod, only: versH2H2cia, strictboundcia, H2orthopara_mixture - 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 pres ! 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 amagat - 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 ios - - character(LEN=*), parameter :: fmat11 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)" - character(LEN=*), parameter :: fmat18 = "(A12,A3,A5,F10.6,F10.4,I7,F7.3,E10.3,F5.3)" - - character*20 bleh - double precision blah, Ttemp, ztemp - integer nres - - integer ind - - ztemp=temp - if ((H2orthopara_mixture .eq. "hot")) then - if (ztemp .gt. 7000.) then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too high for this H2-H2 CIA dataset (>7000K, Hot Jupiter case)' - endif !is_master - call abort_physic("interpolateH2H2", "temperatures are too high for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too high for this H2-H2 CIA dataset (Hot Jupiter case)' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 7000 ' - ! print*,' ... for H2-H2 CIA dataset ... ' - ! print*,' ... we assume we know what you are doing ... ' - ! print*,'*********************************************************' - !endif !is_master - ztemp = 7000. - endif !strictboundcia - endif !(temp .gt. 7000.) - else ! if not "hot" - if(ztemp.gt.400)then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too high for this H2-H2 CIA dataset. If you ' - print*,'really want to run simulations with hydrogen at T > 400 K, contact' - print*,'Robin Wordsworth [rwordsworth@uchicago.edu].' - endif !is_master - call abort_physic("interpolateH2H2", "temperatures are too high for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too high for this H2-H2 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 400 ' - ! print*,' ... for H2-H2 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-H2 CIA dataset. If you ' - print*,'really want to run simulations with hydrogen at T < 40 K, contact' - print*,'Robin Wordsworth [rwordsworth@uchicago.edu].' - endif ! is_master - call abort_physic("interpolateH2H2", "temperatures are too low for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too low for this H2-H2 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 40 ' - ! print*,' ... for H2-H2 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").and. (temp .gt. 3000.)) - amagat = (273.15/temp)*(pres/101325.0) - - if(firstcall)then ! called by sugas_corrk only - if (is_master) print*,'----------------------------------------------------' - if (is_master) print*,'Initialising H2-H2 continuum from HITRAN database...' - - ! 1.1 Open the ASCII files and set nS according to version - ! Only two possible versions for now : 2011 or 2018 (sanity check in inifis_mod) - if (versH2H2cia.eq.2011) then - nS = 2428 - if (H2orthopara_mixture.eq."normal") then - dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia' - else if (H2orthopara_mixture.eq."equilibrium") then - dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_2011.cia' - else if (H2orthopara_mixture.eq."hot") then - dt_file=TRIM(datadir)//'/continuum_data/H2-H2_2011_extended.cia' - ns = 800 - endif - else if (versH2H2cia.eq.2018 .and. H2orthopara_mixture.eq."normal") then - nS = 9600 - dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2018.cia' - else if (versH2H2cia.eq.2018 .and. H2orthopara_mixture.eq."equilibrium") then - nS = 10860 - dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_0-15000cm-1_40-400K.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 interpolateH2H2' - 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_physic("interpolateH2H2", "Unable to read file", 1) - else - - if(versH2H2cia.eq.2011) then - if (is_master) then - write(*,*) '... You are using H2-H2 CIA from 2011 but you should use more recent data available on HITRAN !' - write(*,*) '... (Especially if you are running a giant planet atmosphere)' - write(*,*) '... Just find out the H2-H2 CIA from 2018, put it in your datadir and have a look at interpolateH2H2.F90 ! .' - endif - endif - - do iT=1,nT - - ! Only two possibles values for now : 2011 or 2018 (sanity check in inifis_mod) - if(versH2H2cia.eq.2011) then - read(33,fmat11) bleh,blah,blah,nres,Ttemp - else if (versH2H2cia.eq.2018) then - read(33,fmat18) bleh,bleh,bleh,blah,blah,nres,Ttemp - endif - - 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 interpolateH2H2.F90' - endif - call abort_physic("interpolateH2H2", "Resolution given does not match with nS value", 1) - 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*,'interpolateH2H2: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' pressure ',pres,' Pa' - endif - endif - - - - 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*amagat**2 ! convert to m^-1 - - !print*,'We have ',amagat,' amagats of H2' - !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 interpolateH2H2 Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateH2He.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateH2He.F90 (revision 4054) +++ (revision ) @@ -1,212 +1,0 @@ - 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 - call abort_physic("interpolateH2He", "temperatures are too high for this CIA dataset", 1) - 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 - call abort_physic("interpolateH2He", "temperatures are too high for this CIA dataset", 1) - 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 - call abort_physic("interpolateH2He", "temperatures are too low for this CIA dataset", 1) - 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_physic("interpolateH2He", "Unable to read file", 1) - 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 - call abort_physic("interpolateH2He", "Resolution given does not match with nS value", 1) - 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 Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateH2O_self_foreign.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateH2O_self_foreign.F90 (revision 4054) +++ (revision ) @@ -1,176 +1,0 @@ - subroutine interpolateH2O_self_foreign(wn,temp,presS,presF,abcoef,firstcall,ind) - -!================================================================== -! -! Purpose -! ------- -! Calculates the H2O continuum self and foreign opacity -! using lookup tables from MT_CKD version 3.3 -! -! Authors -! ------- -! M. Turbet (2020) -! Adapted from the file of R. Wordsworth -! -!================================================================== - - 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 presS ! self-pressure (Pascals) - double precision presF ! foreign (air) pressure (Pascals) - - ! output - double precision abcoef ! absorption coefficient (m^-1) - - integer nS,nT, iT - parameter(nS=2001) ! number of wavenumbers - parameter(nT=39) ! number of temperatures - - double precision kB - parameter(kB=1.380649e-23) - - double precision amagatS, amagatF, abcoefS, abcoefF, Nmolec - double precision wn_arr(nS) - double precision temp_arr(nT) - double precision abs_arrS(nS,nT) - double precision abs_arrF(nS,nT) - double precision data_tmp(nT) - - character*43 dummy_var1 - character*46 dummy_var2 - - integer nres - - double precision Ttemp - - character(LEN=*), parameter :: fmat1 = "(A43,I6,F6.1)" - character(LEN=*), parameter :: fmat2 = "(A46,I6,F6.1)" - - - integer k,ind - logical firstcall - - save wn_arr, temp_arr, abs_arrS, abs_arrF !read by master - - character*100 dt_file - integer strlen,ios - - amagatS=(273.15/temp)*(presS/101325.0) - amagatF=(273.15/temp)*(presF/101325.0) - - if(firstcall)then ! called by sugas_corrk only - if (is_master) print*,'----------------------------------------------------' - if (is_master) print*,'Initialising H2O continuum from MT_CKD data...' - -! 1.1 Open the ASCII files - -!$OMP MASTER - - ! self broadening - dt_file=TRIM(datadir)//'/continuum_data/H2O-H2O_continuum_MT_CKD3.3.cia' - open(34,file=dt_file,form='formatted',status='old',iostat=ios) - if (ios.ne.0) then ! file not found - if (is_master) then - write(*,*) 'Error from interpolateH2O_cont SELF' - 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(*,*) ' Check that there is a H2O-H2O_continuum_MT_CKD3.3.cia' - write(*,*)'Continuum file available here:' - write(*,*)' https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/continuum_data/pre-2025_continua/' - endif - call abort_physic("interpolateH2OH2O", "Unable to read file", 1) - else - do iT=1,nT - read(34,fmat1) dummy_var1,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*,'Adjust nS value in interpolateH2O_MTCKD...F90' - endif - call abort_physic("interpolateH2OH2O", "Resolution given does not match with nS value", 1) - endif - temp_arr(iT)=Ttemp - !write(*,*) 'H2O_H2O' - do k=1,nS - read(34,*) wn_arr(k),abs_arrS(k,it) - !write(*,*) nres,temp_arr(iT),wn_arr(k),abs_arrS(k,it) - end do - end do - endif - close(34) - - ! foreign broadening - dt_file=TRIM(datadir)//'/continuum_data/H2O-AIR_continuum_MT_CKD3.3.cia' - open(35,file=dt_file,form='formatted',status='old',iostat=ios) - if (ios.ne.0) then ! file not found - if (is_master) then - write(*,*) 'Error from interpolateH2O_cont FOREIGN' - 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(*,*)' Check that there is a H2O-AIR_continuum_MT_CKD3.3.cia' - write(*,*)'Continuum file available here:' - write(*,*)' https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/continuum_data/pre-2025_continua/' - endif - call abort_physic("interpolateH2OH2O", "Unable to read file", 1) - else - do iT=1,nT - read(35,fmat2) dummy_var2,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*,'Adjust nS value in interpolateH2O_MTCKD...F90' - endif - call abort_physic("interpolateH2OH2O", "Resolution given does not match with nS value", 1) - endif - temp_arr(iT)=Ttemp - !write(*,*) 'H2O_AIR' - do k=1,nS - read(35,*) wn_arr(k),abs_arrF(k,it) - !write(*,*) nres,temp_arr(iT),wn_arr(k),abs_arrF(k,it) - end do - end do - endif - close(35) - if (is_master) then - print*,'interpolateH2O_MTCKDcont: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' H2O pressure ',presS,' Pa' - print*,' air pressure ',presF,' Pa' - endif -!$OMP END MASTER -!$OMP BARRIER - - endif - - call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS,ind) - !print*,'MTCKD - self absorption is ',abcoefS,' cm^2 molecule^-1' - - call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF,ind) - !print*,'MTCKD - foreign absorption is ',abcoefF,' cm^2 molecule^-1' - -! abcoefS and abcoefF are in cm-1 amagat-2 -! First we multiply by amagat^2 - abcoef = abcoefS*amagatS + abcoefF*amagatF - abcoef = abcoef*amagatS -! Then we convert cm-1 in m-1 - abcoef = 100.*abcoef - - !print*,'MTCKD_v3.3 : So the total absorption is ',abcoef,' m^-1' - !print*,'for PH2O/PN2/TEMP/wn=',presS,presF,temp,wn - - return - end subroutine interpolateH2O_self_foreign - - Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateHeCH4.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateHeCH4.F90 (revision 4054) +++ (revision ) @@ -1,176 +1,0 @@ -subroutine interpolateHeCH4(wn,temp,presHe,presCH4,abcoef,firstcall,ind) - - !================================================================== - ! - ! Purpose - ! ------- - ! Calculates the He-CH4 CIA opacity, using a lookup table from - ! HITRAN (2018) - ! - ! Authors - ! ------- - ! R. Wordsworth (2011) - ! - !================================================================== - - use callkeys_mod, only: 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 presCH4 ! CH4 partial pressure (Pascals) - double precision presHe ! He partial pressure (Pascals) - - ! output - double precision abcoef ! absorption coefficient (m^-1) - - integer nS,nT - parameter(nS=1000) - 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 amagatCH4 - 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 !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(ztemp.gt.350)then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too high for this He-CH4 CIA dataset.' - print*,'Please run mixed He-CH4 atmospheres below T = 350 K.' - endif - call abort_physic("interpolateHeCH4", "temperatures are too high for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too high for this He-CH4 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 350 ' - ! print*,' ... for He-CH4 CIA dataset ... ' - ! print*,' ... we assume we know what you are doing ... ' - ! print*,'*********************************************************' - !endif - ztemp = 350 - endif - elseif(ztemp.lt.40)then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too low for this He-CH4 CIA dataset.' - print*,'Please run mixed He-CH4 atmospheres above T = 40 K.' - endif - call abort_physic("interpolateHeCH4", "temperatures are too low for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too low for this He-CH4 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 40 ' - ! print*,' ... for He-CH4 CIA dataset ... ' - ! print*,' ... we assume we know what you are doing ... ' - ! print*,'*********************************************************' - !endif - ztemp = 40 - endif - endif - - amagatCH4 = (273.15/temp)*(presCH4/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 He-CH4 continuum from HITRAN database...' - - ! 1.1 Open the ASCII files - dt_file=TRIM(datadir)//'/continuum_data/CH4-He_2018.cia' - -!$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 interpolateHeCH4' - 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/He-CH4_2018.cia is there.' - endif - call abort_physic("interpolateHeCH4", "Unable to read file", 1) - 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 interpolateHeCH4.F90' - endif - call abort_physic("interpolateHeCH4", "Resolution given does not match with nS value", 1) - 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*,'interpolateHeCH4: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' CH4 partial pressure ',presCH4,' Pa' - print*,' and He partial pressure ',presHe,' Pa' - endif - endif - - 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*amagatCH4*amagatHe ! convert to m^-1 - -! print*,'We have ',amagatCH4,' amagats of CH4' -! 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 interpolateHeCH4 - Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateN2H2.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateN2H2.F90 (revision 4054) +++ (revision ) @@ -1,158 +1,0 @@ -subroutine interpolateN2H2(wn,temp,presN2,presH2,abcoef,firstcall,ind) - - !================================================================== - ! - ! Purpose - ! ------- - ! Calculates the N2-H2 CIA opacity, using a lookup table from - ! HITRAN (2011) - ! - ! Authors - ! ------- - ! R. Wordsworth (2011) - ! - !================================================================== - - use callkeys_mod, only: 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 presN2 ! N2 partial pressure (Pascals) - double precision presH2 ! H2 partial pressure (Pascals) - - ! output - double precision abcoef ! absorption coefficient (m^-1) - - integer nS,nT - parameter(nS=1914) - 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 amagatN2 - double precision amagatH2 - 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 !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(ztemp.gt.400)then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too high for this N2-H2 CIA dataset. If you ' - print*,'really want to run simulations with hydrogen at T > 400 K,' - print*,'find relevant data.' - endif !is_master - call abort_physic("interpolateN2H2", "temperatures are too high for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too high for this N2-H2 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 4000 ' - ! print*,' ... for N2-H2 CIA dataset ... ' - ! print*,' ... we assume we know what you are doing ... ' - ! print*,'*********************************************************' - !endif !is_master - ztemp = 400. - endif !strictboundcia - endif - - amagatN2 = (273.15/temp)*(presN2/101325.0) - amagatH2 = (273.15/temp)*(presH2/101325.0) - - if(firstcall)then ! called by sugas_corrk only - if (is_master) print*,'----------------------------------------------------' - if (is_master) print*,'Initialising N2-H2 continuum from HITRAN database...' - - ! 1.1 Open the ASCII files - dt_file=TRIM(datadir)//'/continuum_data/N2-H2_2011.cia' - -!$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 interpolateN2H2' - 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/N2-H2_2011.cia is there.' - endif - call abort_physic("interpolateN2H2","Unable to read file",1) - 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 interpolateN2H2.F90' - endif - call abort_physic("interpolateN2H2", "Resolution given does not match with nS value", 1) - 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*,'interpolateN2H2: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' N2 partial pressure ',presN2,' Pa' - print*,' and H2 partial pressure ',presH2,' Pa' - endif - endif - - 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*amagatN2*amagatH2 ! convert to m^-1 - -! print*,'We have ',amagatN2,' amagats of N2' -! print*,'and ',amagatH2,' amagats of H2' -! 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 interpolateN2H2 - Index: trunk/LMDZ.GENERIC/libf/phygeneric/interpolateN2N2.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/interpolateN2N2.F90 (revision 4054) +++ (revision ) @@ -1,153 +1,0 @@ -subroutine interpolateN2N2(wn,temp,pres,abcoef,firstcall,ind) - - !================================================================== - ! - ! Purpose - ! ------- - ! Calculates the N2-N2 CIA opacity, using a lookup table from - ! HITRAN (2011) - ! - ! Authors - ! ------- - ! R. Wordsworth (2011) - ! - !================================================================== - - use callkeys_mod, only: 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 pres ! pressure (Pascals) - integer :: ind - - ! output - double precision abcoef ! absorption coefficient (m^-1) - - integer nS,nT - parameter(nS=582) - 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 amagat - 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 !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 - - ztemp=temp - - if(ztemp.gt.400)then - if (strictboundcia) then - if (is_master) then - print*,'Your temperatures are too high for this N2-N2 CIA dataset. If you ' - print*,'really want to run simulations with hydrogen at T > 400 K,' - print*,'find relevant data.' - endif !is_master - call abort_physic("interpolateN2N2", "temperatures are too high for this CIA dataset", 1) - else - !if (is_master) then - ! print*,'Your temperatures are too high for this N2-?2 CIA dataset' - ! print*,'you have chosen strictboundcia = ', strictboundcia - ! print*,'*********************************************************' - ! print*,' we allow model to continue but with temp = 4000 ' - ! print*,' ... for N2-N2 CIA dataset ... ' - ! print*,' ... we assume we know what you are doing ... ' - ! print*,'*********************************************************' - !endif !is_master - ztemp = 400. - endif !strictboundcia - endif - - - amagat = (273.15/temp)*(pres/101325.0) - - if(firstcall)then ! called by sugas_corrk only - if (is_master) print*,'----------------------------------------------------' - if (is_master) print*,'Initialising N2-N2 continuum from HITRAN database...' - - ! 1.1 Open the ASCII files - dt_file=TRIM(datadir)//'/continuum_data/N2-N2_2011.cia' - -!$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 interpolateN2N2' - 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/N2-N2_norm_2011.cia is there.' - endif - call abort_physic("interpolateN2N2", "Unable to read file", 1) - 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 interpolateN2N2.F90' - endif - call abort_physic("interpolateN2N2", "Resolution given does not match with nS value", 1) - 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*,'interpolateN2N2: At wavenumber ',wn,' cm^-1' - print*,' temperature ',temp,' K' - print*,' pressure ',pres,' Pa' - endif - endif - 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*amagat**2 ! convert to m^-1 -! abcoef=0. - -! print*,'We have ',amagat,' amagats of N2' -! 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 interpolateN2N2 - Index: trunk/LMDZ.GENERIC/libf/phygeneric/optci.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/optci.F90 (revision 4054) +++ trunk/LMDZ.GENERIC/libf/phygeneric/optci.F90 (revision 4055) @@ -11,10 +11,9 @@ use radinc_h, only: L_LEVELS, L_NLAYRAD, L_NSPECTI, L_NGAUSS, & L_NLEVRAD, L_REFVAR, naerkind - use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig + use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,glat_ig use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2, & igas_CH4, igas_CO2, igas_O2 use comcstfi_mod, only: g, r, mugaz - use callkeys_mod, only: kastprof,continuum,graybody,varspec, & - generic_continuum_database + use callkeys_mod, only: kastprof,continuum,graybody,varspec use recombin_corrk_mod, only: corrk_recombin, gasi_recomb use tpindex_mod, only: tpindex @@ -98,6 +97,4 @@ integer igas, jgas - - integer interm logical :: firstcall=.true. @@ -126,10 +123,4 @@ endif ! of if (firstcall) - !! AS: to save time in computing continuum (see bilinearbig) - IF (.not.ALLOCATED(indi)) THEN - ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx)) - indi = -9999 ! this initial value means "to be calculated" - ENDIF - !======================================================================= ! Determine the total gas opacity throughout the column, for each @@ -191,5 +182,4 @@ ! include continua if necessary - if(generic_continuum_database)then T_cont = dble(TMID(k)) do igas=1,ngasmx @@ -241,113 +231,4 @@ enddo ! igas=1,ngasmx - - else ! generic_continuum_database - - wn_cont = dble(wnoi(nw)) - T_cont = dble(TMID(k)) - do igas=1,ngasmx - - if(gfrac(igas).eq.-1)then ! variable - p_cont = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol - elseif(varspec) then - p_cont = dble(PMID(k)*scalep*FRACVAR(igas,k)*(1.-QVAR(k))) - else - p_cont = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k))) - endif - - dtemp=0.0d0 - if(igas.eq.igas_N2)then - - interm = indi(nw,igas,igas) - call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) - indi(nw,igas,igas) = interm - - elseif(igas.eq.igas_H2)then - - ! first do self-induced absorption - interm = indi(nw,igas,igas) - call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) - indi(nw,igas,igas) = interm - - ! then cross-interactions with other gases - do jgas=1,ngasmx - if(varspec) then - p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k))) - else - p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k))) - endif - dtempc = 0.0d0 - if(jgas.eq.igas_N2)then - interm = indi(nw,igas,jgas) - call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indi(nw,igas,jgas) = interm - elseif(jgas.eq.igas_CO2)then - interm = indi(nw,igas,jgas) - call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indi(nw,igas,jgas) = interm - elseif(jgas.eq.igas_He)then - interm = indi(nw,igas,jgas) - call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indi(nw,igas,jgas) = interm - endif - dtemp = dtemp + dtempc - enddo - - elseif(igas.eq.igas_CH4)then - - ! first do self-induced absorption - interm = indi(nw,igas,igas) - call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm) - indi(nw,igas,igas) = interm - - ! then cross-interactions with other gases - do jgas=1,ngasmx - if(varspec) then - p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k))) - else - p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k))) - endif - dtempc = 0.0d0 - if(jgas.eq.igas_H2)then - interm = indi(nw,igas,jgas) - call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indi(nw,igas,jgas) = interm - elseif(jgas.eq.igas_CO2)then - interm = indi(nw,igas,jgas) - call interpolateCO2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indi(nw,igas,jgas) = interm - elseif(jgas.eq.igas_He)then - interm = indi(nw,igas,jgas) - call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indi(nw,igas,jgas) = interm - endif - dtemp = dtemp + dtempc - enddo - - elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then - ! Compute self and foreign (with air) continuum of H2O - p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air! - interm = indi(nw,igas,igas) - call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3 - indi(nw,igas,igas) = interm - - endif - - DCONT = DCONT + dtemp - - enddo - - ! Oobleck test - !rho = PMID(k)*scalep / (TMID(k)*286.99) - !if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then - ! DCONT = rho * 0.125 * 4.6e-4 - !elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then - ! DCONT = 1000*dpr(k) * 1.0 * 4.6e-4 / g - ! DCONT = rho * 1.0 * 4.6e-4 - !elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then - ! DCONT = rho * 0.125 * 4.6e-4 - !endif - - endif ! generic_continuum_database DCONT = DCONT*dz(k) Index: trunk/LMDZ.GENERIC/libf/phygeneric/optcv.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/optcv.F90 (revision 4054) +++ trunk/LMDZ.GENERIC/libf/phygeneric/optcv.F90 (revision 4055) @@ -10,10 +10,10 @@ use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND - use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig + use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,glat_ig use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2, & igas_CH4, igas_CO2, igas_O2 use comcstfi_mod, only: g, r, mugaz use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec, & - generic_continuum_database,rayleigh + rayleigh use recombin_corrk_mod, only: corrk_recombin, gasv_recomb use tpindex_mod, only: tpindex @@ -107,6 +107,4 @@ integer igas, jgas - integer interm - logical :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) @@ -120,10 +118,4 @@ firstcall=.false. endif ! of if (firstcall) - - !! AS: to save time in computing continuum (see bilinearbig) - IF (.not.ALLOCATED(indv)) THEN - ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx)) - indv = -9999 ! this initial value means "to be calculated" - ENDIF !======================================================================= @@ -211,5 +203,4 @@ ! include continua if necessary - if(generic_continuum_database)then T_cont = dble(TMID(k)) do igas=1,ngasmx @@ -261,107 +252,4 @@ enddo ! igas=1,ngasmx - - else ! generic_continuum_database - - - wn_cont = dble(wnov(nw)) - T_cont = dble(TMID(k)) - do igas=1,ngasmx - - if(gfrac(igas).eq.-1)then ! variable - p_cont = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol - elseif(varspec) then - p_cont = dble(PMID(k)*scalep*FRACVAR(igas,k)*(1.-QVAR(k))) - else - p_cont = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k))) - endif - - dtemp=0.0 - if(igas.eq.igas_N2)then - - interm = indv(nw,igas,igas) -! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) - indv(nw,igas,igas) = interm - ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible - - elseif(igas.eq.igas_H2)then - - ! first do self-induced absorption - interm = indv(nw,igas,igas) - call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) - indv(nw,igas,igas) = interm - - ! then cross-interactions with other gases - do jgas=1,ngasmx - if(varspec) then - p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k))) - else - p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k))) - endif - dtempc = 0.0 - if(jgas.eq.igas_N2)then - interm = indv(nw,igas,jgas) - call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indv(nw,igas,jgas) = interm - ! should be irrelevant in the visible - elseif(jgas.eq.igas_CO2)then - interm = indv(nw,igas,jgas) - call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indv(nw,igas,jgas) = interm - ! might not be relevant in the visible - elseif(jgas.eq.igas_He)then - interm = indv(nw,igas,jgas) - call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indv(nw,igas,jgas) = interm - endif - dtemp = dtemp + dtempc - enddo - - elseif(igas.eq.igas_CH4)then - - ! first do self-induced absorption - interm = indv(nw,igas,igas) - call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm) - indv(nw,igas,igas) = interm - - ! then cross-interactions with other gases - do jgas=1,ngasmx - if(varspec) then - p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k))) - else - p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k))) - endif - dtempc = 0.0d0 - if(jgas.eq.igas_H2)then - interm = indv(nw,igas,jgas) - call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indv(nw,igas,jgas) = interm - elseif(jgas.eq.igas_CO2)then - interm = indv(nw,igas,jgas) - call interpolateCO2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indv(nw,igas,jgas) = interm - ! might not be relevant in the visible - elseif(jgas.eq.igas_He)then - interm = indv(nw,igas,jgas) - call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) - indv(nw,igas,jgas) = interm - endif - dtemp = dtemp + dtempc - enddo - - elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then - ! Compute self and foreign (with air) continuum of H2O - p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air! - interm = indv(nw,igas,igas) - call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3 - indv(nw,igas,igas) = interm - - endif - - DCONT = DCONT + dtemp - - enddo - - endif ! generic_continuum_database DCONT = DCONT*dz(k) Index: trunk/LMDZ.GENERIC/libf/phygeneric/radcommon_h.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/radcommon_h.F90 (revision 4054) +++ trunk/LMDZ.GENERIC/libf/phygeneric/radcommon_h.F90 (revision 4055) @@ -73,9 +73,4 @@ REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90 !$OMP THREADPRIVATE(blami,blamv) - - !! AS: introduced to avoid doing same computations again for continuum - INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indi - INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indv -!$OMP THREADPRIVATE(indi,indv) !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011 Index: trunk/LMDZ.GENERIC/libf/phygeneric/sugas_corrk.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phygeneric/sugas_corrk.F90 (revision 4054) +++ trunk/LMDZ.GENERIC/libf/phygeneric/sugas_corrk.F90 (revision 4055) @@ -34,5 +34,5 @@ use ioipsl_getin_p_mod, only: getin_p use callkeys_mod, only: varactive,varfixed,graybody,callgasvis,& - continuum,generic_continuum_database + continuum use tracer_h, only : nqtot, moderntracdef, is_recomb, noms use recombin_corrk_mod, only: su_recombin, & @@ -62,6 +62,4 @@ double precision testcont ! for continuum absorption initialisation - - integer :: dummy if (.not. moderntracdef) use_premix=.true. ! Added by JVO for compatibility with 'old' traceur.def @@ -717,5 +715,5 @@ ! Initialise the continuum absorption data if(continuum)then - if(generic_continuum_database)then + do igas=1,ngasmx ! we loop on all pairs of molecules that have data available ! data can be downloaded from https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/continuum_data/ @@ -806,61 +804,4 @@ enddo ! igas=1,ngasmx - else ! generic_continuum_database tag - do igas=1,ngasmx - if (igas .eq. igas_N2) then - - dummy = -9999 - call interpolateN2N2(100.D+0,250.D+0,17500.D+0,testcont,.true.,dummy) - - elseif (igas .eq. igas_H2) then - - ! first do self-induced absorption - dummy = -9999 - call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.,dummy) - ! then cross-interactions with other gases - do jgas=1,ngasmx - if (jgas .eq. igas_N2) then - dummy = -9999 - call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.,dummy) - elseif (jgas .eq. igas_CO2) then - dummy = -9999 - call interpolateCO2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.,dummy) - elseif (jgas .eq. igas_He) then - dummy = -9999 - call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy) - endif - enddo - - - elseif (igas .eq. igas_CH4) then - - ! first do self-induced absorption - dummy = -9999 - call interpolateCH4CH4(600.0,66.7,400.0,testcont,.true.,dummy) - ! then cross-interactions with other gases - do jgas=1,ngasmx - if (jgas .eq. igas_H2) then - dummy = -9999 - call interpolateH2CH4(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy) - elseif (jgas .eq. igas_CO2) then - dummy = -9999 - call interpolateCO2CH4(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.,dummy) - elseif (jgas .eq. igas_He) then - dummy = -9999 - call interpolateHeCH4(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy) - endif - enddo - - - elseif (igas .eq. igas_H2O) then - - ! Compute self and foreign (with air) continuum of H2O - dummy = -9999 - call interpolateH2O_self_foreign(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.,dummy) - - endif - - enddo ! igas=1,ngasmx - endif ! generic_continuum_database tag endif ! continuum flag