Changeset 2245 for trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90

Timestamp:

Feb 24, 2020, 6:55:50 PM (5 years ago)

Author:

jvatant

Message:

+ Add a 'versH2H2cia' int key in callphys that allows two values (2011 or 2018) to

deal with updated HITRAN file (for interpolateH2H2.F90) from 2018 that includes the
H2H2 dimer from Fletcher et al. 2018, useful for giant planets.
Retrocompatibility is ok, default value to 2011.

--JVO

File:

• trunk/LMDZ.GENERIC/libf/phystd/interpolateH2H2.F90 (modified) (5 diffs)

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

@@ -6 +6 @@
 !     -------
 !     Calculates the H2-H2 CIA opacity, using a lookup table from
-!     HITRAN (2011)
+!     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
       use datafile_mod, only: datadir
 
@@ -27 +31 @@
 
       integer nS,nT
-      parameter(nS=2428)
       parameter(nT=10)
 
@@ -36 +39 @@
 
       double precision amagat
-      double precision wn_arr(nS)
       double precision temp_arr(nT)
-      double precision abs_arr(nS,nT)
+      
+      double precision, dimension(:),   allocatable :: wn_arr
+      double precision, dimension(:,:), allocatable :: abs_arr
 
       integer k,iT
       logical firstcall
 
-      save wn_arr, temp_arr, abs_arr !read by master
+      save nS, wn_arr, temp_arr, abs_arr !read by master
 
       character*100 dt_file
-      integer strlen,ios
+      integer ios
 
-      character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)"
+      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
@@ -69 +74 @@
          print*,'Initialising H2-H2 continuum from HITRAN database...'
 
-!     1.1 Open the ASCII files
-         dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia'
+!     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
+           dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia'
+           nS = 2428
+         else if (versH2H2cia.eq.2018) then
+           dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2018.cia'
+           nS = 9600
+         endif
+
+         if(.not.allocated(wn_arr))  allocate(wn_arr(nS)) 
+         if(.not.allocated(abs_arr)) allocate(abs_arr(nS,nT)) 
 
 !$OMP MASTER
@@ -80 +95 @@
            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-H2_norm_2011.cia is there.'
+           write(*,*) 'Also check that the continuum data continuum_data/H2-H2_norm_2011.cia or H2-H2_norm_2018.cia is there.'
            call abort
          else
+         
+         if(versH2H2cia.eq.2011) 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_norm_2018.cia, put it in your datadir and have a look at interpolateH2H2.F90 ! .'
+         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
 
-               read(33,fmat1) bleh,blah,blah,nres,Ttemp
                if(nS.ne.nres)then
                   print*,'Resolution given in file: ',trim(dt_file)