subroutine exponesph(h,ps,cofx,lmax,nc,ncp1,imu) c c*********************************************************************** cccc c subroutine exponesph c c purpose- c 1. calculate the slant path scattering and total optical depths c for a particular line of sight. This is used in the c integration of the output beam. c 2. Then calculate the trapezoidal integration factors using the c calculated slant path optical depths. c c 3. calculate the iteration integration factors for a flat c atmosphere for eq. 4.11 and k functions of eq 4.6. c This is done the same way as in routine expone. c c method- c 1. using the fine layer scattering and total optical depths in c ps and dxs, and the altitude of the layer boundaries, h, c use simple spherical trig. to calc. the ratio of the distance c through a layer in a spherical atmosphere to that in a flat c atmosphere. The result is stored in the arrays ttsf_sph and c ttf_sph. c 2. The integration factors are computed using extmu=layer scattering c optical thickness * transmission to from level to top of atmo. c Integral = 0.5*Source(1)*Transmission(1)* c scattering optical thickness of layer 1 + c Sum from i=2 to nc of 0.5*(Source(i)*Transmission(i) c +Source(i+1)*Transmission(i+1))* c scattering optical thickness of layer i + c 0.5*Source(ncp1)*Transmission(ncp1) c scattering optical thickness of layer ncp1 c 3. If the scan angle is equal to 0.0 degrees, then the quantities c are computed quickly using a flat atmosphere. Again, note c that the integration factors used in the iteration process c are for a flat atmosphere. Only the final integration of the c out-going beam is done in a spherical atmosphere. c c calling sequence- c call exponesph(h,ps,cofx,lmax,nc,ncp1,imu) c c variable type i/o description c -------- ---- --- ----------- c c h(487) r*8 i height from earth center(fraction of r) c ps(487) r*8 i rayleigh opt. thickness of each layer c of the standard atmosphere c cofx(4,487) r*8 i spline interpolation coeff. c lmax i*4 i number of layers in the standard atmosphere c nc i*4 i # of layers from top of c atmosphere to reflecting surface c ncp1 i*4 i nc+1 c imu i*4 i # of scan angles c passed through common blocks c dxs(487) r*8 i total opt. thickness of each layer c of the standard atmosphere c xs(487) r*8 i log of the total opt. depth to a level c of the standard atmosphere c ttl(202) r*8 i log of the total vertical optical depth c emu(15) r*8 i cosine of the scan angle c extmu(202,15) r*8 o trapezoidal integration factor for each level c c internal arrays of note c dtsp_sph(202) r*8 avg. scattering optical thickness of each c spherical layer c dtts_sph(202) r*8 avg. scattering optical thickness of each c spherical half layer at mid points between c layers c transph(202) r*8 transmission from level to the top of the c spherical atmosphere c c external references c dexpk1 (subroutine) c omerf (function) c c common areas referenced c consts c eks c emm c es c prints c thkns c chpmn c csphout c kmat c crefdir c depolt c c last modified 03/14/95...dave flittner c purpose: set pressure scale height used in gravity correction c to rayleigh scattering od. Create new variable pscaleforg and c pass in common block consts. c*********************************************************************** cccc implicit none c input integer*4 lmax,nc,ncp1,imu real*8 h(487),ps(487),cofx(4,487) c internal integer*4 i,j,lmaxm1,im,im1,im2,jk,k,ip1,ip2,ii,iii real*8 sintheta,pc,pcpc,dist1,dist2,ddist,ratio,dum1,dum2,dum3 real*8 ttsf_sph(487),ttf_sph(487) real*8 tts_sph(202) real*8 dtts_sph(202),dtsp_sph(202),transph(202) real*8 chpp,chxx,thta,sum,thetap real*8 ekary(202,5),extmux(202,15) c c common block area c include "chpmn.inc" include "thkns.inc" include "csphout.inc" include "prints.inc" include "emm.inc" include "kmat.inc" include "eks.inc" include "consts.inc" include "es.inc" include "crefdir.inc" c c new statement include "depolt.inc" c end of new statement c c external functions used real*8 omerf external omerf c end of portion lifted from expone ! def c c loop for each scan angle c do 199 im=1,imu c now loop using the simple trig. relations c calc. the tangent height in terms of earth radii using the scan angle c at the surface sintheta=(1.0d0-emu(im)*emu(im)) if(sintheta.le.0.0d0)then c do not need to correct, for sphericity, so will use the flat atmo. case do i=1,ncp1 extmu(i,im)=dtsp(i)*dexp(-tt(i)) enddo refdir(im)=dexp(-tt(ncp1)) goto 199 else if(sintheta.gt.1.0d0)then sintheta=1.0d0 pc=h(lmax) else sintheta=dsqrt(sintheta) pc=h(lmax)*sintheta endif c start at top of atmosphere and work way down along the line of sight c use chapman function to do the first layer which ranges from h(1) to c infinity. Note that the approximation used to evaluate the Chapman c function for dum1=1.0 (scan angle of 0.0 deg.) returns a value less c than that in a flat atmosphere. dum1=pc/h(1) dum1=dsqrt(1.0d0-dum1**2) chpp=omerf(sqchp*dum1,dum1,chpn) chxx=omerf(sqchx*dum1,dum1,chxn) ttsf_sph(1)=ps(1)*chpp ttf_sph(1)=(dxs(1)-ps(1))*chxx+ps(1)*chpp c now start the loop to create the fine arrays pcpc=pc*pc dist1=dsqrt(h(1)*h(1)-pcpc) do i=2,lmax dist2=dsqrt(h(i)*h(i)-pcpc) ddist=dist1-dist2 ratio=ddist/(h(i-1)-h(i)) ttsf_sph(i)=ttsf_sph(i-1)+ps(i)*ratio ttf_sph(i)=ttf_sph(i-1)+dxs(i)*ratio dist1=dist2 enddo sum=0.0d0 c now spline log of slant total od. as a function of the log of veritcal od. do i=1,lmax ttf_sph(i)=log(ttf_sph(i)) enddo call splset(xs,ttf_sph,cofx,lmax) c c*****for each layer of model atmosphere find two layers in the c*****standard atmosphere straddling it and find total slant optical c*****path to reach each layer of model atmosphere- c*****using spline interpolation and then compute the slant tranmission c*****and store in transph c transph(1) = 1.0 j = 2 lmaxm1 = lmax - 1 do 110 i = 1, lmaxm1 if (ttl(j) .le. xs(i) .or. ttl(j) .gt. xs(i+1)) go to 110 100 dum1 = xs(i+1) - ttl(j) dum2 = ttl(j) - xs(i) dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1 dum2**2 + cofx(4,i)) dum3 = dexp(dum3) transph(j) = dexp(-dum3) j = j + 1 if (j .gt. ncp1) go to 120 if (ttl(j) .gt. xs(i) .and. ttl(j) .le. xs(i+1)) go to 100 110 continue c c*****find values for the bottom layer by extrapolation if it was not c*****obtained during the above interpolation c 120 if (ttl(ncp1) .le. xs(lmax)) go to 125 dum1 = xs(lmax) - ttl(ncp1) dum2 = ttl(ncp1) - xs(lmaxm1) i = lmaxm1 dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1 dum2**2 + cofx(4,i)) dum3 = dexp(dum3) transph(ncp1) = dexp(-dum3) 125 continue refdir(im)=transph(ncp1) c c now need to calculate the spherical scattering optical thickness of c each layer. Already have the scattering optical depth along the slant path, c so can use a spline of log(tts_sph) versus xs and interpolate to the c points ttl and use the same method to compute the average scattering c optical thickness of each layer as is done in dtaus do i=1,lmax ttsf_sph(i)=log(ttsf_sph(i)) enddo call splset(xs,ttsf_sph,cofx,lmax) c now interpolate to ttl to get slant scattering od. c c*****for each layer of model atmosphere find two layers in the c*****standard atmosphere straddling it and find slant scattering optical c*****path to reach each layer of model atmosphere- c*****using spline interpolation and store in tts_sph c tts_sph(1) = 0.0 j = 2 lmaxm1 = lmax - 1 do 140 i = 1, lmaxm1 if (ttl(j) .le. xs(i) .or. ttl(j) .gt. xs(i+1)) go to 140 130 dum1 = xs(i+1) - ttl(j) dum2 = ttl(j) - xs(i) dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1 dum2**2 + cofx(4,i)) dum3 = dexp(dum3) tts_sph(j) = dum3 j = j + 1 if (j .gt. ncp1) go to 150 if (ttl(j) .gt. xs(i) .and. ttl(j) .le. xs(i+1)) go to 130 140 continue c c*****find values for the bottom layer by extrapolation if it was not c*****obtained during the above interpolation c 150 if (ttl(ncp1) .le. xs(lmax)) go to 155 dum1 = xs(lmax) - ttl(ncp1) dum2 = ttl(ncp1) - xs(lmaxm1) i = lmaxm1 dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1 dum2**2 + cofx(4,i)) dum3 = dexp(dum3) tts_sph(ncp1) = dum3 155 continue c c now calc. dtts_sph and dtsp_sph c dtts_sph(1) = 0.5*tts_sph(2) dtsp_sph(1) = 0.5*tts_sph(2) do i=2,nc dtsp_sph(i) = 0.5*(tts_sph(i+1) - tts_sph(i-1)) dtts_sph(i) = 0.5*(tts_sph(i) - tts_sph(i-1)) enddo dtts_sph(ncp1) = 0.5*(tts_sph(ncp1)-tts_sph(nc)) dtsp_sph(ncp1) = dtts_sph(ncp1) c c now calc. the integration factors c dum1=emu(im) do i=1,ncp1 extmu(i,im)=dtsp_sph(i)*transph(i)*dum1 enddo 199 continue c c portion lifted from expone !def c do 400 i=1,ncp1 c initialize ekary do 200 j=1,5 ekary(i,j)=0.d0 200 continue c c calculate w(t)dt*attenuation factor (see eqn. 3.6,6.7) c do 205 j=1,imu extmux(i,j)=extmu(i,j)/emu(j) !def 205 continue if (i.eq.1) go to 255 im1=i-1 c do 210 j=1,5 ekary(im1,j)=0. 210 continue c if(i.eq.2) go to 235 im2=i-2 c c calculates the elements of k1 matrix (see eqn 4.11) c evaluate k1(1),k1(2) of eqn 4.7 for current value of c argument (tt(i)-tt(j)) c do 225 j=1,im2 call dexpk(tt(i),tt(j)) c c ekary array contains k matrix*w(t)dt c do 220 jk=1,5 ekary(j,jk)=dtsp(j)*eek(jk) 220 continue if(i.lt.ncp1) go to 225 c new statements if(ipol.eq.0)then ek4(j)=(e(2)+e(4))*dtsp(j) ek5(j)=sq2*(e(2)-e(4))*dtsp(j) else ek4(j)=((1.0d0+2.0d0*q)*e(2)+e(4))*dtsp(j) ek5(j)=delp*(e(2)-e(4))*dtsp(j) endif c end of new statements 225 continue c call dexpk(tt(i),tt(im1)) do 230 k=1,5 ekary(im1,k)=dtts(im1)*eek(k) 230 continue c if(i.lt.ncp1) go to 235 c new statements if(ipol.eq.0)then ek4(nc)=dtsp(nc)*(e(2)+e(4)) ek5(nc)=sq2*dtsp(nc)*(e(2)-e(4)) else ek4(nc)=((1.0d0+2.0d0*q)*e(2)+e(4))*dtsp(nc) ek5(nc)=delp*(e(2)-e(4))*dtsp(nc) endif c end of new statements 235 call dexpk1(tt(i),tt(im1)) c do 240 k=1,5 ekary(im1,k)=ekary(im1,k)+dtts(i)*eek(k) 240 continue c do 250 k=1,5 ekary(i,k)=dtts(i)*eek(k) 250 continue if(i.lt.ncp1) go to 255 c new statemnets if(ipol.eq.0)then ek4(ncp1)=dtsp(ncp1)*(odan(2)+odan(4)) ek5(ncp1)=sq2*dtsp(ncp1)*(odan(2)-odan(4)) else ek4(ncp1)=((1.0d0+2.0d0*q)*odan(2)+odan(4))*dtsp(ncp1) ek5(ncp1)=delp*(odan(2)-odan(4))*dtsp(ncp1) endif c end of new statements go to 300 255 ip1=i+1 call dexpk1(tt(i),tt(ip1)) c do 260 k=1,5 ekary(i,k)=ekary(i,k)+dtts(ip1)*eek(k) 260 continue c do 270 k=1,5 ekary(ip1,k)=dtts(ip1)*eek(k) 270 continue c if(i.eq.nc) go to 300 ip2=i+2 call dexpk(tt(i),tt(ip1)) c do 280 k=1,5 ekary(ip1,k)=ekary(ip1,k)+dtts(ip2)*eek(k) 280 continue c do 295 j=ip2,ncp1 call dexpk(tt(i),tt(j)) c do 290 k=1,5 ekary(j,k)=dtsp(j)*eek(k) 290 continue 295 continue c c 300 continue do 320 ii = 1,202 do 310 iii = 1,5 nek(ii,iii,i) = ekary(ii,iii) 310 continue 320 continue 400 continue c c optionally print calculated exponential terms c (zstar matrix * w(t)dt see eqn 6.3) c if (jprint(3).ne.0) then write (33,6120) (j,ek4(j),ek5(j) ,j=1,ncp1) do 495 i=1,imu thta=dacos(emu(i))/cnvrt write(33,7120) thta,(j,extmux(j,i),j=1,ncp1,10) write(33,7130) ncp1,extmux(ncp1,i) 495 continue endif return 6120 format(1h1,t50,'debug print out for expone',///, 1 1h ,2(' j',4x,'ek4(j)',6x,'ek5(j)',6x),//,2(1x,i3,2d12.4,3x)) 7120 format(1x,//,1x,'theta=',f6.2,' *** attenuation factors/*** ',//, 1 1x,'level',5x,'extmu/mu',/,(1x,i5,3x,d12.6)) 7130 format(1x,i5,3x,d12.6) c end of portion lifted from expone.f !def end