source: trunk/LMDZ.GENERIC/libf/phystd/interpolateH2Ocont_CKD.F90 @ 848

Last change on this file since 848 was 716, checked in by rwordsworth, 12 years ago

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

Note to Martians: I've changed the default values of 'sedimentation' and
'co2cond' in inifis to false. Both these are defined in the standard deftank
callphys.def file, so there shouldn't be any compatibility problems.
File size: 8.2 KB
Line 
1     subroutine interpolateH2Ocont_CKD(wn,temp,presS,presF,abcoef,firstcall)
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
43      logical firstcall
44
45      save wn_arr, temp_arr, abs_arrS, abs_arrF
46
47      character*100 dt_file
48      integer strlen,ios
49
50      amagatS=(273.15/temp)*(presS/101325.0)
51      amagatF=(273.15/temp)*(presF/101325.0)
52
53      if(firstcall)then ! called by sugas_corrk only
54         print*,'----------------------------------------------------'
55         print*,'Initialising H2O continuum from MT_CKD data...'
56
57!     1.1 Open the ASCII files
58
59         ! nu array
60         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_NU.dat'
61         open(33,file=dt_file,form='formatted',status='old',iostat=ios)
62         if (ios.ne.0) then        ! file not found
63           write(*,*) 'Error from interpolateH2O_cont'
64           write(*,*) 'Data file ',trim(dt_file),' not found.'
65           write(*,*)'Check that your path to datagcm:',trim(datadir)
66           write(*,*)' is correct. You can change it in callphys.def with:'
67           write(*,*)' datadir = /absolute/path/to/datagcm'
68           write(*,*)' Also check that there is a continuum_data/H2O_CONT_NU.dat there.'
69           call abort
70         else
71            do k=1,nS
72               read(33,*) wn_arr(k)
73            enddo
74         endif
75         close(33)
76
77         ! self broadening
78         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_SELF.dat'
79         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
80         if (ios.ne.0) then        ! file not found
81           write(*,*) 'Error from interpolateH2O_cont'
82           write(*,*) 'Data file ',trim(dt_file),' not found.'
83           write(*,*)'Check that your path to datagcm:',trim(datadir)
84           write(*,*)' is correct. You can change it in callphys.def with:'
85           write(*,*)' datadir = /absolute/path/to/datagcm'
86           write(*,*)' Also check that there is a continuum_data/H2O_CONT_SELF.dat there.'
87           call abort
88         else
89            do k=1,nS
90               read(34,*) data_tmp
91               abs_arrS(k,1:nT)=data_tmp(1:nT)
92            end do
93         endif
94         close(34)
95
96         ! foreign (N2+O2+Ar) broadening
97         dt_file=TRIM(datadir)//'/continuum_data/H2O_CONT_FOREIGN.dat'
98         open(35,file=dt_file,form='formatted',status='old',iostat=ios)
99         if (ios.ne.0) then        ! file not found
100           write(*,*) 'Error from interpolateH2O_cont'
101           write(*,*) 'Data file ',trim(dt_file),' not found.'
102           write(*,*)'Check that your path to datagcm:',trim(datadir)
103           write(*,*)' is correct. You can change it in callphys.def with:'
104           write(*,*)' datadir = /absolute/path/to/datagcm'
105           write(*,*)' Also check that there is a continuum_data/H2O_CONT_FOREIGN.dat there.'
106           call abort
107         else
108            do k=1,nS
109               read(35,*) data_tmp
110               abs_arrF(k,1:nT)=data_tmp(1:nT)
111            end do
112         endif
113         close(35)
114
115         temp_arr(1)  = 200.
116         temp_arr(2)  = 250.
117         temp_arr(3)  = 300.
118         temp_arr(4)  = 350.
119         temp_arr(5)  = 400.
120         temp_arr(6)  = 450.
121         temp_arr(7)  = 500.
122         temp_arr(8)  = 550.
123         temp_arr(9)  = 600.
124         temp_arr(10) = 650.
125         temp_arr(11) = 700.
126
127         print*,'interpolateH2Ocont: At wavenumber ',wn,' cm^-1'
128         print*,'   temperature ',temp,' K'
129         print*,'   H2O pressure ',presS,' Pa'
130         print*,'   air pressure ',presF,' Pa'
131
132         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
133         print*,'the self absorption is ',abcoefS,' cm^2 molecule^-1'
134
135         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
136         print*,'the foreign absorption is ',abcoefF,' cm^2 molecule^-1'
137
138         print*,'We have ',amagatS,' amagats of H2O vapour'
139         print*,'and ',amagatF,' amagats of air'
140
141         abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
142         abcoef = abcoef*(presS/(presF+presS))      ! take H2O mixing ratio into account
143                                                    ! abs coeffs are given per molecule of H2O
144
145         Nmolec = (presS+presF)/(kB*temp)           ! assume ideal gas
146         print*,'Total number of molecules per m^3 is',Nmolec
147
148         abcoef = abcoef*Nmolec/(100.0**2)          ! convert to m^-1
149         print*,'So the total absorption is ',abcoef,' m^-1'
150         print*,'And optical depth / km : ',1000.0*abcoef
151
152      else
153
154         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
155         call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
156
157         abcoef = abcoefS*amagatS + abcoefF*amagatF
158         abcoef = abcoef*(presS/(presF+presS))
159
160         Nmolec = (presS+presF)/(kB*temp)
161         abcoef = abcoef*Nmolec/(100.0**2)
162
163         if(wn.gt.500 .and. wn.lt.1400)then
164         elseif(wn.gt.2100 .and. wn.lt.3000)then
165         else
166            abcoef = 0.0
167         endif
168
169         ! unlike for Rayleigh scattering, we do not currently weight by the BB function
170         ! however our bands are normally thin, so this is no big deal.
171
172      endif
173
174      return
175    end subroutine interpolateH2Ocont_CKD
176
177
178!-------------------------------------------------------------------------
179      subroutine bilinearH2Ocont(x_arr,y_arr,f2d_arr,x_in,y_in,f)
180!     Necessary for interpolation of continuum data
181
182      implicit none
183
184      integer nX,nY,i,j,a,b
185      parameter(nX=1001)
186      parameter(nY=11)
187
188      real*8 x_in,y_in,x,y,x1,x2,y1,y2
189      real*8 f,f11,f12,f21,f22,fA,fB
190      real*8 x_arr(nX)
191      real*8 y_arr(nY)
192      real*8 f2d_arr(nX,nY)
193     
194      integer strlen
195      character*100 label
196      label='subroutine bilinear'
197
198      x=x_in
199      y=y_in
200
201!     1st check we're within the wavenumber range
202      if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
203         f=0.0D+0
204         return
205      else
206         
207!     in the x (wavenumber) direction 1st
208         i=1
209 10      if (i.lt.(nX+1)) then
210            if (x_arr(i).gt.x) then
211               x1=x_arr(i-1)
212               x2=x_arr(i)
213               a=i-1
214               i=9999
215            endif
216            i=i+1
217            goto 10
218         endif
219      endif
220     
221      if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
222         write(*,*) 'Warning from bilinearH2Ocont:'
223         write(*,*) 'Outside continuum temperature range!'
224         if(y.lt.y_arr(1))then
225            y=y_arr(1)+0.01
226         endif
227         if(y.gt.y_arr(nY))then
228            y=y_arr(nY)-0.01
229         endif
230      else
231
232!     in the y (temperature) direction 2nd
233         j=1
234 20      if (j.lt.(nY+1)) then
235            if (y_arr(j).gt.y) then
236               y1=y_arr(j-1)
237               y2=y_arr(j)
238               b=j-1
239               j=9999
240            endif
241            j=j+1
242            goto 20
243         endif
244      endif
245     
246      f11=f2d_arr(a,b)
247      f21=f2d_arr(a+1,b)
248      f12=f2d_arr(a,b+1)
249      f22=f2d_arr(a+1,b+1)
250     
251!     1st in x-direction
252      fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1)
253      fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1)
254     
255!     then in y-direction
256      f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
257
258      return
259    end subroutine bilinearH2Ocont
Note: See TracBrowser for help on using the repository browser.