-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcalconwater_locg.F
178 lines (154 loc) · 6.99 KB
/
calconwater_locg.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
c http://www.mathworks.com/matlabcentral/fileexchange/25934-fortran-95-interface-to-matlab-api-with-extras
c change *.f to *.F
c change integer m_in,n_in to mwSize m_in,n_in
c change integer nlhs,nrhs to integer*2 nlhs,nrhs
c change mxCreateFull to mxCreateDoubleMatrix
c include fintrf.h
#include "fintrf.h"
subroutine mexFunction(nlhs,plhs,nrhs,prhs)
include '/home/sergio/SPECTRA/FORTRANFILES/max.inc'
mwpointer plhs(*),prhs(*)
mwpointer mxGetPr,mxCreateDoubleMatrix
integer nlhs,nrhs
mwsize mxGetM,mxGetN,mx(13),nx(13)
mwSize m_in,n_in
mwsize o_in,p_in
mwPointer ip,np,fp,fsp,nlap,tp,pp,ppp,ap,wckdp,wlp,zp
mwPointer selfp,forp
real*8 raT(kMaxLayer),raP(kMaxLayer),raPP(kMaxLayer)
real*8 raA(kMaxLayer),raZ(MaxLen),whichlayer
real*8 idgas,nfreq,raFreq(MaxLen),fstep,nlay,ckd
real*8 self,for
c check for proper number of arguments
c SUBROUTINE CALCON23( CON, IDGAS, NFREQ, FREQ, FSTEP, NLAY, T, P,
c $ PARTP, AMNT, ckd, selfmult,formult, WHICHLAYER)
c the call is con = CALCON**( IDGAS, NFREQ, FREQ, FSTEP, NLAY, T, P,
c $ PARTP, AMNT, WHICHCKD, selfmult,formult, WHICHLAYER)
c freq= raFreq= array of wavenumbers = Freq(*)
c T,P,PartP,amnt = arrays of layer variables <= 100 == T(*),P(*) etc
c con = raa(kMaxPts) == CON(*)
c WHICHCKD = integer = 00,21,23
c WHICHLAYER = which gas amt etc to use
c selfmult = parameter to multiply self part with
c formult = parameter to multiply for part with
if (nrhs .ne. 13) then
call mexErrMsgTxt('13 input args required')
endif
if (nlhs .ne. 1) then
call mexErrMsgTxt('1 output arg required')
endif
do ii = 1,nrhs
o_in=mxGetM(prhs(ii))
p_in=mxGetN(prhs(ii))
mx(ii) = o_in
nx(ii) = p_in
end do
c want to check sizes of input array "freq"
m_in=mxGetM(prhs(3))
n_in=mxGetN(prhs(3))
if ((m_in .gt. MaxLen) .or. (n_in .gt. MaxLen)) then
call mexErrMsgTxt('wavevector has to be smaller than MaxLen')
endif
if ((m_in .ne. 1) .and. (n_in .ne. 1)) then
call mexErrMsgTxt('input wavevector needs to be 1d array')
endif
c want to check sizes of input array "T,P,PP,A"
o_in=mxGetM(prhs(6))
p_in=mxGetN(prhs(6))
if ((o_in .gt. kMaxLayer) .or. (p_in .gt. kMaxLayer)) then
call mexErrMsgTxt('p,pp,t,a has to be smaller than kMaxLayer')
endif
if ((o_in .ne. 1) .and. (p_in .ne. 1)) then
call mexErrMsgTxt('p,pp,t,a wavevector needs to be 1d array')
endif
ip = mxGetPr(prhs(1))
np = mxGetPr(prhs(2))
fp = mxGetPr(prhs(3))
fsp = mxGetPr(prhs(4))
nlap = mxGetPr(prhs(5))
tp = mxGetPr(prhs(6))
pp = mxGetPr(prhs(7))
ppp = mxGetPr(prhs(8))
ap = mxGetPr(prhs(9))
wckdp = mxGetPr(prhs(10))
selfp = mxGetPr(prhs(11))
forp = mxGetPr(prhs(12))
wlp = mxGetPr(prhs(13))
c copy right hand arguments to local arrays or variables
c z = calcon**()
c note that in reality nbox and zlen are integers
call mxCopyPtrToReal8(ip, idgas, 1)
call mxCopyPtrToReal8(np, nfreq, 1)
call mxCopyPtrToReal8(fp, raFreq, int(max(mx(3),nx(3))))
call mxCopyPtrToReal8(fsp, fstep, 1)
call mxCopyPtrToReal8(nlap, nlay, 1)
call mxCopyPtrToReal8(tp, raT, int(max(o_in,p_in)))
call mxCopyPtrToReal8(pp, raP, int(max(o_in,p_in)))
call mxCopyPtrToReal8(ppp, raPP, int(max(o_in,p_in)))
call mxCopyPtrToReal8(ap, raA, int(max(o_in,p_in)))
call mxCopyPtrToReal8(wckdp, ckd ,1)
call mxCopyPtrToReal8(selfp, self ,1)
call mxCopyPtrToReal8(forp, for ,1)
call mxCopyPtrToReal8(wlp, whichlayer,1)
c create a matrix for return argument and assign pointers to the
c output parameters z = boxint3(y,nbox,zlen)
plhs(1) = mxCreateDoubleMatrix(m_in,n_in,0)
zp = mxGetPr(plhs(1))
c do the actual computations in a subroutine
c these first few are original and.or modifications of MT_CKD 1
if (int(ckd) .eq. 01) then
print *,'the MT_CKD version = ',int(ckd)
call calcon_mtckd_01_loc(raZ,int(idgas),int(nfreq),raFreq,
$ fstep,int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
elseif (int(ckd) .eq. 02) then
print *,'the MT_CKD_Hannon version = ',int(ckd)
call calcon_mtckd_02_loc(raZ,int(idgas),int(nfreq),raFreq,
$ fstep,int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
elseif (int(ckd) .eq. 03) then
print *,'the MT_CKD_Hannon version = ',int(ckd)
call calcon_mtckd_03_loc(raZ,int(idgas),int(nfreq),raFreq,
$ fstep,int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
elseif (int(ckd) .eq. 04) then
print *,'the MT_CKD_Hannon version = ',int(ckd)
call calcon_mtckd_04_loc(raZ,int(idgas),int(nfreq),raFreq,
$ fstep,int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
elseif (int(ckd) .eq. 05) then
print *,'the MT_CKD_Hannon version = ',int(ckd)
call mexErrMsgTxt('for CKD5 = Tobin, use /home/sergio/tobin5.m')
c call calcon_mtckd_05_loc(raZ,int(idgas),int(nfreq),raFreq,
c $ fstep,int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
elseif (int(ckd) .eq. 06) then
print *,'the MT_CKD_Hannon version = ',int(ckd)
call calcon_mtckd_06_loc(raZ,int(idgas),int(nfreq),raFreq,
$ fstep,int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
c these are MT-CKD 2.5 or mods of that, DONE SEPARATE
c elseif (int(ckd) .eq. 25) then
c print *,'the MT_CKD version = ',int(ckd)
c call calcon_mtckd_25_loc(raZ,int(idgas),int(nfreq),raFreq,
c $ fstep,int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
c these last few are the old CKD models
elseif (int(ckd) .eq. 00) then
print *,'the CKD version = ',int(ckd)
call calcon00_loc(raZ,int(idgas),int(nfreq),raFreq,fstep,
$ int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
elseif (int(ckd) .eq. 21) then
print *,'the CKD version = ',int(ckd)
call calcon21_loc(raZ,int(idgas),int(nfreq),raFreq,fstep,
$ int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
elseif (int(ckd) .eq. 23) then
print *,'the CKD version = ',int(ckd)
call calcon23_loc(raZ,int(idgas),int(nfreq),raFreq,fstep,
$ int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
elseif (int(ckd) .eq. 24) then
print *,'the CKD version = ',int(ckd)
call calcon24_loc(raZ,int(idgas),int(nfreq),raFreq,fstep,
$ int(nlay),raT,raP,raPP,raA,self,for,int(whichlayer))
else
print *,'the CKD version = ',int(ckd)
print *,'code needs <0>,<21>,<23>,<24> or <1>,2,3,4,5 or 6'
call mexErrMsgTxt('please retry')
end if
c copy output which is stored in local array to matrix output
call mxCopyReal8ToPtr(raZ, zp, max(m_in,n_in))
return
end