source: trunk/LMDZ.VENUS/libf/phyvenus/interpolateH2Ocont_CKD.F90 @ 3556

Last change on this file since 3556 was 2560, checked in by slebonnois, 3 years ago

SL: Implementation of SW computation based on generic model. Switch between this new SW module or old module that reads R. Haus tables implemented with a key (solarchoice)

File size: 6.1 KB
Line 
1     subroutine interpolateH2Ocont_CKD(wn,temp,presS,presF,abcoef,firstcall,ind)
2
3!==================================================================
4!     
5!     Purpose
6!     -------
7!     Calculates the H2O continuum opacity, using a lookup table from
8!     Clough (2005).
9!
10!     Authors
11!     -------
12!     R. Wordsworth (2011)
13!     
14!==================================================================
15
16      use datafile_mod, only: datadir
17      implicit none
18
19      ! input
20      double precision wn                 ! wavenumber             (cm^-1)
21      double precision temp               ! temperature            (Kelvin)
22      double precision presS              ! self-pressure          (Pascals)
23      double precision presF              ! foreign (air) pressure (Pascals)
24
25      ! output
26      double precision abcoef             ! absorption coefficient (m^-1)
27
28      integer nS,nT
29      parameter(nS=1001)
30      parameter(nT=11)
31
32      double precision kB
33      parameter(kB=1.3806488e-23)
34
35      double precision amagatS, amagatF, abcoefS, abcoefF, Nmolec
36      double precision wn_arr(nS)
37      double precision temp_arr(nT)
38      double precision abs_arrS(nS,nT)
39      double precision abs_arrF(nS,nT)
40      double precision data_tmp(nT)
41
42      integer k,ind
43      logical firstcall
44
45      save wn_arr, temp_arr, abs_arrS, abs_arrF !read by master
46
47      character*100 dt_file
48      integer strlen,ios
49
50      character(len=100) :: message
51      character(len=10),parameter :: subname="H2OcontCKD"
52
53      amagatS=(273.15/temp)*(presS/101325.0)
54      amagatF=(273.15/temp)*(presF/101325.0)
55
56      if(firstcall)then ! called by sugas_corrk only
57         print*,'----------------------------------------------------'
58         print*,'Initialising H2O continuum from MT_CKD data...'
59
60!     1.1 Open the ASCII files
61
62!$OMP MASTER
63         ! nu array
64         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_NU.dat'
65         open(33,file=dt_file,form='formatted',status='old',iostat=ios)
66         if (ios.ne.0) then        ! file not found
67           write(*,*) 'Error from interpolateH2O_cont'
68           write(*,*) 'Data file ',trim(dt_file),' not found.'
69           write(*,*)'Check that your path to datagcm:',trim(datadir)
70           write(*,*)' is correct. You can change it in callphys.def with:'
71           write(*,*)' datadir = /absolute/path/to/datagcm'
72           message=' Also check that there is a continuum_data/H2O_CONT_NU.dat there.'
73           call abort_physic(subname,message,1)
74         else
75            do k=1,nS
76               read(33,*) wn_arr(k)
77            enddo
78         endif
79         close(33)
80
81         ! self broadening
82         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_SELF.dat'
83         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
84         if (ios.ne.0) then        ! file not found
85           write(*,*) 'Error from interpolateH2O_cont'
86           write(*,*) 'Data file ',trim(dt_file),' not found.'
87           write(*,*)'Check that your path to datagcm:',trim(datadir)
88           write(*,*)' is correct. You can change it in callphys.def with:'
89           write(*,*)' datadir = /absolute/path/to/datagcm'
90           message=' Also check that there is a continuum_data/H2O_CONT_SELF.dat there.'
91           call abort_physic(subname,message,1)
92         else
93            do k=1,nS
94               read(34,*) data_tmp
95               abs_arrS(k,1:nT)=data_tmp(1:nT)
96            end do
97         endif
98         close(34)
99
100         ! foreign (N2+O2+Ar) broadening
101         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_FOREIGN.dat'
102         open(35,file=dt_file,form='formatted',status='old',iostat=ios)
103         if (ios.ne.0) then        ! file not found
104           write(*,*) 'Error from interpolateH2O_cont'
105           write(*,*) 'Data file ',trim(dt_file),' not found.'
106           write(*,*)'Check that your path to datagcm:',trim(datadir)
107           write(*,*)' is correct. You can change it in callphys.def with:'
108           write(*,*)' datadir = /absolute/path/to/datagcm'
109           message=' Also check that there is a continuum_data/H2O_CONT_FOREIGN.dat there.'
110           call abort_physic(subname,message,1)
111         else
112            do k=1,nS
113               read(35,*) data_tmp
114               abs_arrF(k,1:nT)=data_tmp(1:nT)
115            end do
116         endif
117         close(35)
118
119         temp_arr(1)  = 200.
120         temp_arr(2)  = 250.
121         temp_arr(3)  = 300.
122         temp_arr(4)  = 350.
123         temp_arr(5)  = 400.
124         temp_arr(6)  = 450.
125         temp_arr(7)  = 500.
126         temp_arr(8)  = 550.
127         temp_arr(9)  = 600.
128         temp_arr(10) = 650.
129         temp_arr(11) = 700.
130
131         print*,'interpolateH2Ocont: At wavenumber ',wn,' cm^-1'
132         print*,'   temperature ',temp,' K'
133         print*,'   H2O pressure ',presS,' Pa'
134         print*,'   air pressure ',presF,' Pa'
135!$OMP END MASTER
136!$OMP BARRIER
137         
138      endif
139
140      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS,ind)
141!      print*,'the self absorption is ',abcoefS,' cm^2 molecule^-1'
142
143      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF,ind)
144!      print*,'the foreign absorption is ',abcoefF,' cm^2 molecule^-1'
145
146!      print*,'We have ',amagatS,' amagats of H2O vapour'
147!      print*,'and ',amagatF,' amagats of air'
148
149      abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
150      abcoef = abcoef*(presS/(presF+presS))      ! take H2O mixing ratio into account
151                                                    ! abs coeffs are given per molecule of H2O
152
153      Nmolec = (presS+presF)/(kB*temp)           ! assume ideal gas
154!      print*,'Total number of molecules per m^3 is',Nmolec
155
156      abcoef = abcoef*Nmolec/(100.0**2)          ! convert to m^-1
157!      print*,'So the total absorption is ',abcoef,' m^-1'
158!      print*,'And optical depth / km : ',1000.0*abcoef
159
160
161      if(wn.gt.500 .and. wn.lt.1400)then
162      elseif(wn.gt.2100 .and. wn.lt.3000)then
163      else
164         abcoef = 0.0
165      endif
166
167      ! unlike for Rayleigh scattering, we do not currently weight by the BB function
168      ! however our bands are normally thin, so this is no big deal.
169
170
171      return
172    end subroutine interpolateH2Ocont_CKD
173
Note: See TracBrowser for help on using the repository browser.