Changeset 2655


Ignore:
Timestamp:
Mar 30, 2022, 12:19:20 PM (3 years ago)
Author:
gmilcareck
Message:

Major changes to CIA interpolation:
1) Add contribution from CH4 (H2-CH4,He-CH4,CH4-CH4) ;
2) Add some tests before interpolation for H2, He and CH4 ;
3) Add the possibility to choose between a normal or equilibrium ortho:para
fraction for CIA H2;
4) Change "strictboundH2H2cia" to the generic "strictboundcia" for H2,He,CH4.
It can be added for others CIA (N2,H,CO2...) if you want.

Location:
trunk/LMDZ.GENERIC
Files:
3 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r2625 r2655  
    17151715Minor bug fix for aerave_new.F when input wavelenght in data file are in
    17161716descending order.
     1717
     1718== 30/03/2022 == GM
     1719Major changes to CIA interpolation:
     17201) Add contribution from CH4 (H2-CH4,He-CH4,CH4-CH4) ;
     17212) Add some tests before interpolation for H2, He and CH4 ;
     17223) Add the possibility to choose between a normal or equilibrium ortho:para
     1723fraction for CIA H2;
     17244) Change "strictboundH2H2cia" to the generic "strictboundcia" for H2,He,CH4.
     1725It can be added for others CIA (N2,H,CO2...) if you want.

  • trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90

    r2635 r2655  
    1414      logical,save :: strictboundcorrk                                     
    1515!$OMP THREADPRIVATE(strictboundcorrk)
    16       logical,save :: strictboundH2H2cia                                     
    17 !$OMP THREADPRIVATE(strictboundH2H2cia)
     16      logical,save :: strictboundcia                                     
     17!$OMP THREADPRIVATE(strictboundcia)
    1818
    1919      logical,save :: enertest
    … …  
    7676      integer,save :: startype
    7777      integer,save :: versH2H2cia
     78      character(64),save :: H2orthopara_mixture
    7879      integer,save :: nlayaero
    79 !$OMP THREADPRIVATE(iddist,iaervar,iradia,startype,versH2H2cia,nlayaero)
     80!$OMP THREADPRIVATE(iddist,iaervar,iradia,startype,versH2H2cia,H2orthopara_mixture,nlayaero)
    8081      integer,dimension(:),allocatable,save :: aeronlay_choice
    8182!$OMP THREADPRIVATE(aeronlay_choice)

  • trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90

    r2635 r2655  
    259259
    260260     if (is_master) write(*,*) trim(rname)//&
    261        ": prohibit calculations outside H2H2 CIA T grid?"
    262      strictboundH2H2cia=.true. ! default value
    263      call getin_p("strictboundH2H2cia",strictboundH2H2cia)
    264      if (is_master) write(*,*) trim(rname)//&
    265        ": strictboundH2H2cia = ",strictboundH2H2cia
     261       ": prohibit calculations outside CIA T grid?"
     262     strictboundcia=.true. ! default value
     263     call getin_p("strictboundcia",strictboundcia)
     264     if (is_master) write(*,*) trim(rname)//&
     265       ": strictboundcia = ",strictboundcia
    266266
    267267     if (is_master) write(*,*) trim(rname)//&
    … …  
    304304     if (versH2H2cia.ne.2011 .and. versH2H2cia.ne.2018) then
    305305        call abort_physic(rname,"Error: Choose a correct value (2011 or 2018) for versH2H2cia !",1)
    306      endif
     306     endif
     307     
     308     if (is_master) write(*,*) trim(rname)//&
     309       ": CIA - normal or equilibrium ortho-para mixture for H2?"
     310     H2orthopara_mixture = 'normal'
     311     call getin_p("H2orthopara_mixture",H2orthopara_mixture)
     312     if (is_master) write(*,*)trim(rname)//&
     313       ": H2orthopara_mixture = ",trim(H2orthopara_mixture)
    307314
    308315     if (is_master) write(*,*) trim(rname)//&

  • trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90

    r2633 r2655  
     1
    12     subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,ind)
    23
    … …  
    1718!==================================================================
    1819
    19       use callkeys_mod, only: versH2H2cia, strictboundH2H2cia
     20      use callkeys_mod, only: versH2H2cia, strictboundcia, H2orthopara_mixture
    2021      use datafile_mod, only: datadir
    2122
    … …  
    6061
    6162      integer ind
    62 
     63     
    6364      if(temp.gt.400)then
    64          if (strictboundH2H2cia) then
    65             print*,'Your temperatures are too high for this H2-H2 CIA dataset. If you '
    66             print*,'really want to run simulations with hydrogen at T > 400 K, contact'
    67             print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
    68             stop
    69          else
    70             print*,'Your temperatures are too high for this H2-H2 CIA dataset'
    71             print*,'you have chosen strictboundH2H2cia = ', strictboundH2H2cia
    72             print*,'*********************************************************'
    73             print*,' we allow model to continue but with temp = 400          '
    74             print*,'  ...       for H2-H2 CIA dataset          ...           '
    75             print*,'  ... we assume we know what you are doing ...           '
    76             print*,'*********************************************************'
    77             temp = 400
    78          endif
     65        if (strictboundcia) then
     66          print*,'Your temperatures are too high for this H2-H2 CIA dataset. If you '
     67          print*,'really want to run simulations with hydrogen at T > 400 K, contact'
     68          print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
     69          stop
     70        else
     71          print*,'Your temperatures are too high for this H2-H2 CIA dataset'
     72          print*,'you have chosen strictboundcia = ', strictboundcia
     73          print*,'*********************************************************'
     74          print*,' we allow model to continue but with temp = 400          '
     75          print*,'  ...       for H2-H2 CIA dataset          ...           '
     76          print*,'  ... we assume we know what you are doing ...           '
     77          print*,'*********************************************************'
     78          temp = 400
     79        endif
     80      elseif(temp.lt.40)then     
     81        if (strictboundcia) then
     82          print*,'Your temperatures are too low for this H2-H2 CIA dataset. If you '
     83          print*,'really want to run simulations with hydrogen at T < 40 K, contact'
     84          print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
     85          stop
     86        else
     87          print*,'Your temperatures are too low for this H2-H2 CIA dataset'
     88          print*,'you have chosen strictboundcia = ', strictboundcia
     89          print*,'*********************************************************'
     90          print*,' we allow model to continue but with temp = 40           '
     91          print*,'  ...       for H2-H2 CIA dataset          ...           '
     92          print*,'  ... we assume we know what you are doing ...           '
     93          print*,'*********************************************************'
     94          temp = 40
     95        endif           
    7996      endif
    8097
    … …  
    88105         ! Only two possible versions for now : 2011 or 2018 (sanity check in inifis_mod)
    89106         if (versH2H2cia.eq.2011) then
    90            dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia'
    91107           nS = 2428
     108           if (H2orthopara_mixture.eq."normal") then
     109             dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia'
     110           else if (H2orthopara_mixture.eq."equilibrium") then
     111             dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_2011.cia'
     112           endif
    92113         else if (versH2H2cia.eq.2018) then
    93            dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2018.cia'
    94114           nS = 9600
     115           if (H2orthopara_mixture.eq."normal") then
     116             dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2018.cia'
     117           else if (H2orthopara_mixture.eq."equilibrium") then
     118             dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_2018.cia'
     119           endif
    95120         endif
    96121
    … …  
    106131           write(*,*) 'is correct. You can change it in callphys.def with:'
    107132           write(*,*) 'datadir = /absolute/path/to/datagcm'
    108            write(*,*) 'Also check that the continuum data continuum_data/H2-H2_norm_2011.cia or H2-H2_norm_2018.cia is there.'
     133           write(*,*) 'Also check that the continuum data is there.'
    109134           call abort
    110135         else
    … …  
    113138           write(*,*) '... You are using H2-H2 CIA from 2011 but you should use more recent data available on HITRAN !'
    114139           write(*,*) '... (Especially if you are running a giant planet atmosphere)'
    115            write(*,*) '... Just find out the H2-H2_norm_2018.cia, put it in your datadir and have a look at interpolateH2H2.F90 ! .'
     140           write(*,*) '... Just find out the H2-H2 CIA from 2018, put it in your datadir and have a look at interpolateH2H2.F90 ! .'
    116141         endif
    117142
    … …  
    149174
    150175      endif
    151 
     176     
     177           
     178           
    152179         call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
    153180

  • trunk/LMDZ.GENERIC/libf/phystd/interpolateH2He.F90

    r1315 r2655  
    1414!==================================================================
    1515
     16      use callkeys_mod, only: H2orthopara_mixture, strictboundcia
    1617      use datafile_mod, only: datadir
    1718
    … …  
    5758
    5859      integer ind
    59  
     60      
    6061      if(temp.gt.400)then
    61          print*,'Your temperatures are too high for this H2-He CIA dataset.'
    62          print*,'Please run mixed H2-He atmospheres below T = 400 K.'     
    63          stop
     62        if (strictboundcia) then
     63          print*,'Your temperatures are too high for this H2-He CIA dataset.'
     64          print*,'Please run mixed H2-He atmospheres below T = 400 K.'     
     65          stop
     66        else
     67          print*,'Your temperatures are too high for this H2-He CIA dataset'
     68          print*,'you have chosen strictboundcia = ', strictboundcia
     69          print*,'*********************************************************'
     70          print*,' we allow model to continue but with temp = 400          '
     71          print*,'  ...       for H2-He CIA dataset          ...           '
     72          print*,'  ... we assume we know what you are doing ...           '
     73          print*,'*********************************************************'
     74          temp = 400
     75        endif
     76      elseif(temp.lt.40)then
     77        if (strictboundcia) then
     78          print*,'Your temperatures are too low for this H2-He CIA dataset.'
     79          print*,'Please run mixed H2-He atmospheres above T = 40 K.'     
     80          stop
     81        else
     82          print*,'Your temperatures are too low for this H2-He CIA dataset'
     83          print*,'you have chosen strictboundcia = ', strictboundcia
     84          print*,'*********************************************************'
     85          print*,' we allow model to continue but with temp = 40           '
     86          print*,'  ...       for H2-He CIA dataset          ...           '
     87          print*,'  ... we assume we know what you are doing ...           '
     88          print*,'*********************************************************'
     89          temp = 40
     90        endif
    6491      endif
    6592
    … …  
    7299
    73100!     1.1 Open the ASCII files
    74          dt_file=TRIM(datadir)//'/continuum_data/H2-He_norm_2011.cia'
     101         if (H2orthopara_mixture.eq."normal") then
     102           dt_file=TRIM(datadir)//'/continuum_data/H2-He_norm_2011.cia'
     103         else if (H2orthopara_mixture.eq."equilibrium") then
     104           dt_file=TRIM(datadir)//'/continuum_data/H2-He_eq_2011.cia'
     105         endif
    75106         
    76107!$OMP MASTER
    … …  
    82113           write(*,*) 'is correct. You can change it in callphys.def with:'
    83114           write(*,*) 'datadir = /absolute/path/to/datagcm'
    84            write(*,*) 'Also check that the continuum data continuum_data/H2-He_norm_2011.cia is there.'
     115           write(*,*) 'Also check that the continuum data is there.'
    85116           call abort
    86117         else

  • trunk/LMDZ.GENERIC/libf/phystd/optci.F90

    r2582 r2655  
    1212                      L_NLEVRAD, L_REFVAR, naerkind
    1313  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
    14   use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2
     14  use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2, igas_CH4
    1515  use comcstfi_mod, only: g, r, mugaz
    1616  use callkeys_mod, only: kastprof,continuum,graybody
    … …  
    208208                 enddo
    209209
     210              elseif(igas.eq.igas_CH4)then
     211
     212                 ! first do self-induced absorption
     213                 interm = indi(nw,igas,igas)
     214                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     215                 indi(nw,igas,igas) = interm
     216
     217                 ! then cross-interactions with other gases
     218                 do jgas=1,ngasmx
     219                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
     220                    dtempc  = 0.0d0
     221                    if(jgas.eq.igas_H2)then
     222                       interm = indi(nw,igas,jgas)
     223                       call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     224                       indi(nw,igas,jgas) = interm
     225                    elseif(jgas.eq.igas_He)then
     226                       interm = indi(nw,igas,jgas)
     227                       call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     228                       indi(nw,igas,jgas) = interm
     229                    endif
     230                    dtemp = dtemp + dtempc
     231                 enddo
     232
    210233              elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
    211234                 ! Compute self and foreign (with air) continuum of H2O

  • trunk/LMDZ.GENERIC/libf/phystd/optcv.F90

    r2582 r2655  
    1111  use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND
    1212  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
    13   use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2
     13  use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2, igas_CH4
    1414  use comcstfi_mod, only: g, r, mugaz
    1515  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis
    … …  
    215215                    dtemp = dtemp + dtempc
    216216                 enddo
     217                 
     218              elseif(igas.eq.igas_CH4)then
     219
     220                 ! first do self-induced absorption
     221                 interm = indv(nw,igas,igas)
     222                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     223                 indv(nw,igas,igas) = interm
     224
     225                 ! then cross-interactions with other gases
     226                 do jgas=1,ngasmx
     227                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
     228                    dtempc  = 0.0d0
     229                    if(jgas.eq.igas_H2)then
     230                       interm = indv(nw,igas,jgas)
     231                       call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     232                       indv(nw,igas,jgas) = interm
     233                    elseif(jgas.eq.igas_He)then
     234                       interm = indv(nw,igas,jgas)
     235                       call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     236                       indv(nw,igas,jgas) = interm
     237                    endif
     238                    dtemp = dtemp + dtempc
     239                 enddo
    217240
    218241              elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then

  • trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90

    r2585 r2655  
    3030      use comcstfi_mod, only: mugaz
    3131      use gases_h, only: gnom, ngasmx, &
    32                          igas_H2, igas_H2O, igas_He, igas_N2
     32                         igas_H2, igas_H2O, igas_He, igas_N2, igas_CH4
    3333      use ioipsl_getin_p_mod, only: getin_p
    3434      use callkeys_mod, only: varactive,varfixed,graybody,callgasvis,&
    … …  
    737737               endif
    738738            enddo
     739               
     740         
     741         elseif (igas .eq. igas_CH4) then
     742         
     743            ! first do self-induced absorption
     744            dummy = -9999
     745            call interpolateCH4CH4(600.0,66.7,400.0,testcont,.true.,dummy)
     746            ! then cross-interactions with other gases
     747            do jgas=1,ngasmx
     748               if (jgas .eq. igas_H2) then
     749                  dummy = -9999
     750                  call interpolateH2CH4(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
     751               elseif (jgas .eq. igas_He) then
     752                  dummy = -9999
     753                  call interpolateHeCH4(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
     754               endif
     755            enddo
     756           
    739757
    740758         elseif (igas .eq. igas_H2O) then
Note: See TracChangeset for help on using the changeset viewer.