!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%%
***s/p integ2
*
subroutine integ2(y,r,con,s,a,b,c,n,my,mr,nk,base) 8
*
#include "impnone.cdk"
integer n,my,mr,nk
real y(my,nk),r(mr,nk),con,s(n,nk),a(n,nk),b(n,nk),c(n,nk)
logical base
*
*author
* r. benoit (august 93) adapted from s/r intodif
*
*revision
* 001 G. Pellerin (May 03) - IBM Conversion
* - calls to vslog routine (from massvp4 library)
*
*object
* to resolve differential equation of the 1st order
* d y / d s = con * s**alpha*r with boundary condition at
* s(nk) if base=.true. or at s(1) if base=.false.
*
*arguments
*
* - output -
* y result
*
* - input -
* r right-hand-side (mr,nk)
* con constant muliplier
* alpha exponent of pre-factor of the differential equation
* s sigma levels
* a work space (n,nk)
* b work space (n,nk)
* c work space (n,nk)
* n horizontal dimension
* my 1st dimension of y
* mr 1st dimension of r
* nk vertical dimension
* base .true. for boundary condition at s(nk)
* .false. for boundary condition at s(1)
*
*notes
* y(*,1) or y(*,nk) must be initialized
* based accordingly as y and r cannot share the same space.
*
**
*
real q1,q2,q3,xp,xo,xm,ex,aa,bb,cc,dd,det,ak,bk,ck,co2
real temp(n)
integer j,k,l
*
* calcul des a,b,c
*
*
* cas 1: k = 1
*
k = 1
do j=1,n
xp=s(j,2)
xm=s(j,1)
temp(j)=xp/xm
enddo
call vslog(temp,temp,n)
do j=1,n
xp=s(j,2)
xm=s(j,1)
xo=(xp+xm)*0.5
*gpp q1=alog(xp/xm)
q1=temp(j)
q2=(xp-xm)
q3=(xp*xp-xm*xm)*0.5
q3=q3-xo*(2.0*q2-xo*q1)
q2=q2-xo*q1
aa=xm-xo
bb=xp-xo
cc=aa*aa
dd=bb*bb
det=aa*dd-bb*cc
a(j,k)=(dd*q2-bb*q3)/det*0.5
c(j,k)=(aa*q3-cc*q2)/det*0.5
b(j,k)=q1*0.5-a(j,k)-c(j,k)
aa=a(j,1)
bb=b(j,1)*0.25
cc=c(j,1)
b(j,1)=aa+bb*(1.0+(s(j,3)-s(j,2))/(s(j,3)-s(j,1)))
c(j,1)=cc+bb*(1.0+(s(j,3)-s(j,1))/(s(j,3)-s(j,2)))
a(j,1)=-bb*((s(j,2)-s(j,1))*(s(j,2)-s(j,1)))/
$ ((s(j,3)-s(j,2))*(s(j,3)-s(j,1)))
enddo
*
* cas 2: k = nk
*
k = nk
do j=1,n
xp=s(j,nk)
xm=s(j,nk-1)
temp(j)=xp/xm
enddo
call vslog(temp,temp,n)
do j=1,n
xp=s(j,nk)
xm=s(j,nk-1)
xo=(xp+xm)*0.5
*gpp q1=alog(xp/xm)
q1=temp(j)
q2=(xp-xm)
q3=(xp*xp-xm*xm)*0.5
q3=q3-xo*(2.0*q2-xo*q1)
q2=q2-xo*q1
aa=xm-xo
bb=xp-xo
cc=aa*aa
dd=bb*bb
det=aa*dd-bb*cc
a(j,k)=(dd*q2-bb*q3)/det*0.5
c(j,k)=(aa*q3-cc*q2)/det*0.5
b(j,k)=q1*0.5-a(j,k)-c(j,k)
aa=a(j,nk)
bb=b(j,nk)*0.25
cc=c(j,nk)
b(j,nk)=cc+bb*(1.0+(s(j,nk-1)-s(j,nk-2))/
$ (s(j,nk)-s(j,nk-2)))
a(j,nk)=aa+bb*(1.0+(s(j,nk)-s(j,nk-2))/
$ (s(j,nk-1)-s(j,nk-2)))
c(j,nk)=-bb*((s(j,nk)-s(j,nk-1))*(s(j,nk)-s(j,nk-1)))/
$ ((s(j,nk-1)-s(j,nk-2))*(s(j,nk)-s(j,nk-2)))
enddo
*
* cas 1: k > 1 and < nk
*
do k=2,nk-1
do j=1,n
xp=s(j,k+1)
xm=s(j,k-1)
temp(j)=xp/xm
enddo
call vslog(temp,temp,n)
do j=1,n
xp=s(j,k+1)
xo=s(j,k)
xm=s(j,k-1)
*gpp q1=alog(xp/xm)
q1=temp(j)
q2=(xp-xm)
q3=(xp*xp-xm*xm)*0.5
q3=q3-xo*(2.0*q2-xo*q1)
q2=q2-xo*q1
aa=xm-xo
bb=xp-xo
cc=aa*aa
dd=bb*bb
det=aa*dd-bb*cc
a(j,k)=(dd*q2-bb*q3)/det*0.5
c(j,k)=(aa*q3-cc*q2)/det*0.5
b(j,k)=q1*0.5-a(j,k)-c(j,k)
enddo
enddo
*
*
* integration
*
if (base) then
*
* y(nk) est initialise
*
do 2 j=1,n
ak=-2.0*con*a(j,nk)
bk=-2.0*con*b(j,nk)
ck=-2.0*con*c(j,nk)
2 y(j,nk-1)=ak*r(j,nk-1)+bk*r(j,nk)+ck*r(j,nk-2)+y(j,nk)
do 3 k=nk-2,1,-1
do 3 j=1,n
ak=-2.0*con*a(j,k+1)
bk=-2.0*con*b(j,k+1)
ck=-2.0*con*c(j,k+1)
3 y(j,k)=ak*r(j,k)+bk*r(j,k+1)+ck*r(j,k+2)+y(j,k+2)
else
*
* y(1) est initialise
*
do 4 j=1,n
ak=2.0*con*a(j,1)
bk=2.0*con*b(j,1)
ck=2.0*con*c(j,1)
4 y(j,2)=bk*r(j,1)+ck*r(j,2)+ak*r(j,3)+y(j,1)
do 5 k=3,nk,1
do 5 j=1,n
ak=2.0*con*a(j,k-1)
bk=2.0*con*b(j,k-1)
ck=2.0*con*c(j,k-1)
5 y(j,k)=ak*r(j,k-2)+bk*r(j,k-1)+ck*r(j,k)+y(j,k-2)
endif
*
*
return
end