source: trunk/LMDZ.GENERIC/libf/phystd/interpolateH2O_self_foreign.F90 @ 3523

Last change on this file since 3523 was 2662, checked in by dbardet, 3 years ago

Using of is_master keyword to allow nly the master processor to write the messages and avoid too heavy output files
File size: 6.2 KB
Line 
1     subroutine interpolateH2O_self_foreign(wn,temp,presS,presF,abcoef,firstcall,ind)
2
3!==================================================================
4!     
5!     Purpose
6!     -------
7!     Calculates the H2O continuum self and foreign opacity
8!     using lookup tables from MT_CKD version 3.3
9!
10!     Authors
11!     -------
12!     M. Turbet (2020)
13!     Adapted from the file of R. Wordsworth
14!     
15!==================================================================
16
17      use datafile_mod, only: datadir
18      use mod_phys_lmdz_para, only : is_master
19
20      implicit none
21
22      ! input
23      double precision wn                 ! wavenumber             (cm^-1)
24      double precision temp               ! temperature            (Kelvin)
25      double precision presS              ! self-pressure          (Pascals)
26      double precision presF              ! foreign (air) pressure (Pascals)
27
28      ! output
29      double precision abcoef             ! absorption coefficient (m^-1)
30
31      integer nS,nT, iT
32      parameter(nS=2001)    ! number of wavenumbers
33      parameter(nT=39)      ! number of temperatures
34
35      double precision kB
36      parameter(kB=1.3806488e-23)
37
38      double precision amagatS, amagatF, abcoefS, abcoefF, Nmolec
39      double precision wn_arr(nS)
40      double precision temp_arr(nT)
41      double precision abs_arrS(nS,nT)
42      double precision abs_arrF(nS,nT)
43      double precision data_tmp(nT)
44     
45      character*43 dummy_var1
46      character*46 dummy_var2
47
48      integer nres
49
50      double precision Ttemp
51
52      character(LEN=*), parameter :: fmat1 = "(A43,I6,F6.1)"
53      character(LEN=*), parameter :: fmat2 = "(A46,I6,F6.1)"
54
55
56      integer k,ind
57      logical firstcall
58
59      save wn_arr, temp_arr, abs_arrS, abs_arrF !read by master
60
61      character*100 dt_file
62      integer strlen,ios
63
64      amagatS=(273.15/temp)*(presS/101325.0)
65      amagatF=(273.15/temp)*(presF/101325.0)
66
67      if(firstcall)then ! called by sugas_corrk only
68         if (is_master) print*,'----------------------------------------------------'
69         if (is_master) print*,'Initialising H2O continuum from MT_CKD data...'
70
71!     1.1 Open the ASCII files
72
73!$OMP MASTER
74
75         ! self broadening
76         dt_file=TRIM(datadir)//'/continuum_data/H2O-H2O_continuum_MT_CKD3.3.cia'
77         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
78         if (ios.ne.0) then        ! file not found
79           if (is_master) then
80             write(*,*) 'Error from interpolateH2O_cont SELF'
81             write(*,*) 'Data file ',trim(dt_file),' not found.'
82             write(*,*) 'Check that your path to datagcm:',trim(datadir)
83             write(*,*) ' is correct. You can change it in callphys.def with:'
84             write(*,*) ' datadir = /absolute/path/to/datagcm'
85             write(*,*) ' Check that there is a H2O-H2O_continuum_MT_CKD3.3.cia'
86             write(*,*)'Continuum file available here:'
87             write(*,*)' https://www.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/'
88           endif
89           call abort
90         else
91           do iT=1,nT
92             read(34,fmat1) dummy_var1,nres,Ttemp
93             if(nS.ne.nres)then
94               if (is_master) then
95                 print*,'Resolution given in file: ',trim(dt_file)
96                 print*,'is ',nres,', which does not match nS.'
97                 print*,'Adjust nS value in interpolateH2O_MTCKD...F90'
98               endif
99               stop
100             endif
101             temp_arr(iT)=Ttemp
102             !write(*,*) 'H2O_H2O'
103             do k=1,nS
104               read(34,*) wn_arr(k),abs_arrS(k,it)
105               !write(*,*) nres,temp_arr(iT),wn_arr(k),abs_arrS(k,it)
106             end do
107           end do
108         endif
109         close(34)
110           
111         ! foreign broadening
112         dt_file=TRIM(datadir)//'/continuum_data/H2O-AIR_continuum_MT_CKD3.3.cia'
113         open(35,file=dt_file,form='formatted',status='old',iostat=ios)
114         if (ios.ne.0) then        ! file not found
115           if (is_master) then
116             write(*,*) 'Error from interpolateH2O_cont FOREIGN'
117             write(*,*) 'Data file ',trim(dt_file),' not found.'
118             write(*,*)'Check that your path to datagcm:',trim(datadir)
119             write(*,*)' is correct. You can change it in callphys.def with:'
120             write(*,*)' datadir = /absolute/path/to/datagcm'
121             write(*,*)' Check that there is a H2O-AIR_continuum_MT_CKD3.3.cia'
122             write(*,*)'Continuum file available here:'
123             write(*,*)' https://www.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/'
124           endif
125           call abort
126         else
127           do iT=1,nT
128             read(35,fmat2) dummy_var2,nres,Ttemp
129             if(nS.ne.nres)then
130               if (is_master) then
131                 print*,'Resolution given in file: ',trim(dt_file)
132                 print*,'is ',nres,', which does not match nS.'
133                 print*,'Adjust nS value in interpolateH2O_MTCKD...F90'
134               endif
135               stop
136             endif
137             temp_arr(iT)=Ttemp
138             !write(*,*) 'H2O_AIR'
139             do k=1,nS
140               read(35,*) wn_arr(k),abs_arrF(k,it)
141               !write(*,*) nres,temp_arr(iT),wn_arr(k),abs_arrF(k,it)
142             end do
143           end do
144         endif
145         close(35)
146         if (is_master) then
147           print*,'interpolateH2O_MTCKDcont: At wavenumber ',wn,' cm^-1'
148           print*,'   temperature ',temp,' K'
149           print*,'   H2O pressure ',presS,' Pa'
150           print*,'   air pressure ',presF,' Pa'
151         endif
152!$OMP END MASTER
153!$OMP BARRIER
154         
155      endif
156
157      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS,ind)
158      !print*,'MTCKD - self absorption is ',abcoefS,' cm^2 molecule^-1'
159
160      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF,ind)
161      !print*,'MTCKD - foreign absorption is ',abcoefF,' cm^2 molecule^-1'
162
163! abcoefS and abcoefF are in cm-1 amagat-2
164! First we multiply by amagat^2
165      abcoef = abcoefS*amagatS + abcoefF*amagatF
166      abcoef = abcoef*amagatS
167! Then we convert cm-1 in m-1
168      abcoef = 100.*abcoef
169     
170      !print*,'MTCKD_v3.3 : So the total absorption is ',abcoef,' m^-1'
171      !print*,'for PH2O/PN2/TEMP/wn=',presS,presF,temp,wn
172
173      return
174    end subroutine interpolateH2O_self_foreign
175
176
Note: See TracBrowser for help on using the repository browser.