c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/firstzr.f,v 2.22 2000/08/30 16:38:09 flittner Exp $ subroutine firstzr(h,ps,z,thenot,fracin,lmax,ncp1) c c*********************************************************************** c cccc c subroutine firstzr c c purpose- c correct the optical depths calculated for the layers of a plane c parallel atmosphere for the sphericity of the actual atmosphere c c method- c 1.using the optical depths of the layers of the standard atmos. c calculate the slant optical path of the solar rays thru them c using the chapman function. Also account for refraction of rays. c 2.then using spline interpolation calculate the slant optical c paths in the model atmosphere. c c calling sequence- c call firstzr(h,ps,z,thenot,fracin,fs,lmax,ncp1) c c variables type i/o description c --------- ---- --- ----------- c h(487) r*8 i height of layers of standard atmosphere c from earth center (units of r) c ps(487) r*8 i rayleigh optical thickness of c each layer of std. atmosphere c z(202) r8 o total optical path of solar ray to c reach each layer of the model atmos. c thenot r*8 solar zenith angle c fracin r*8 i layer*r (radius of earth in units of c layer thickness) c fs r*8 o cos(thenot)*z(202) c lmax i*4 i # layers of standard atmos. c ncp1 i*4 i # layers in model atmos. c c external references c omerf c splset 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. cccc c********************************************************************** c c ****************************************************************** c c using the slant path optical thicknesses and chapman function c this subroutine calculates dave's z matrices for primary c scattering. c c ****************************************************************** implicit integer*4(i-n),real*8 (a-h,o-z) include 'parameter.inc' c integer lmax,ncp1 real *8 h(487),ps(487),zs(487),cofx(4,487),fracin,thenot,z(202) include 'consts.cmn' include 'out.cmn' include 'thkns.cmn' include 'chpmn.cmn' include 'energy.cmn' c include 'atm_101.cmn' c local integer maxnz parameter (maxnz=1000) integer nplpts,layindx(2*maxnz),ishadow,iapp,k,mode real*8 losdist(2*maxnz),rpl(2*maxnz+1),mupl(2*maxnz+1) real*8 x(maxnz*2+1),alt,st1,p1,mupt,myzs,y21,y2n,myraycon real*8 thtap(487),thtsfc logical debug,lprt data debug/.false./,lprt/.false./ save debug,lprt c c c*****calculate trig functions to be used in do loops c amuo = dcos(thenot*cnvrt) if (thenot .eq. 90.0) amuo = 0. amuosq = amuo**2 sn = dsqrt(1.-amuosq) c c*****calculate slant opt. path of solar rays to reach c*****layer(i) of standard atmosphere c c c*****find the slant opt. path to reach layer 1 from the top of c*****the atmosphere- using chapman function approximation c i=1 chpp = omerf(sqchp*amuo, amuo, chpn) chxx = omerf(sqchx*amuo, amuo, chxn) zs(i) = chpp*ps(1) + chxx*(dxs(1) - ps(1)) thtap(i)=thenot*dr c do 200 i = 2, lmax c unrefracted ray constant raycon = sn*h(i) c for this altitude, find the mupl(rpl) for the path through the atmosphere iapp=1 tht0=thenot*dr alt=(h(i)-1.)*r call refracray(alt,tht0,r,rfndx,alt_101,n_101, &thtap(i),iapp,nplpts,layindx,losdist,rpl,mupl,x,ishadow,debug) y21=1.0e30 y2n=1.0e30 mode=0 call spline(rpl,mupl,nplpts,y21,y2n,y2_101,mode) myzs=0.0d0 if(lprt)write(6,*)' alt mupt path myzs st1' do j = 1, i alt=h(j)*r call splint(rpl,mupl,y2_101,nplpts,alt,mupt,mode) mupt=max(mupt,-1.0d0) mupt=min(mupt,1.0d0) st1=dsqrt(1.0d0-mupt*mupt) myraycon=h(j)*st1 if(j.eq.1)then c c*****find the slant opt. path to reach layer 1 from the top of c*****the atmosphere- using chapman function approximation amu = myraycon/h(1) amu = dsqrt(1. - amu**2) chpp = omerf(sqchp*amu, amu, chpn) chxx = omerf(sqchx*amu, amu, chxn) myzs = chpp*ps(1) + chxx*(dxs(1) - ps(1)) else path=fracin*(dsqrt((h(j-1)+myraycon)*(h(j-1)-myraycon))- &h(j)*abs(mupt)) myzs=myzs+path*dxs(j) endif if(lprt)then write(6,'(a3,1p6e13.5)')'spl',(h(j)-1.)*r,mupt,path, &myzs,asin(st1)/dr,raycon/myraycon-1. endif enddo zs(i)=myzs 200 continue c c*****fit a spline between xs=log(vert. path) and zs=log(slant path) c*****through the standard atmosphere c do 215 i = 1, lmax 215 zs(i) = dlog(zs(i)) call splset(xs,zs,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 (z) to reach each layer of model atmosphere- c***** using spline interpolation c z(1) = 1.0 j = 2 lmaxm1 = lmax - 1 do 230 i = 1, lmaxm1 if (ttl(j) .le. xs(i) .or. ttl(j) .gt. xs(i+1)) go to 230 220 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)* 1dum2**2 + cofx(4,i)) dum3 = dexp(dum3) z(j) = dexp(-dum3) j = j + 1 if (j .gt. ncp1) go to 240 if (ttl(j) .gt. xs(i) .and. ttl(j) .le. xs(i+1)) go to 220 230 continue c c*****find z for the bottom layer by extrapolation if it was not c*****obtained during the above interpolation c 240 if (ttl(ncp1) .le. xs(lmax)) go to 245 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)* 1dum2**2 + cofx(4,i)) dum3 = dexp(dum3) z(ncp1) = dexp(-dum3) 245 continue c c for this altitude, find the mupl(rpl) for the path through the atmosphere iapp=1 tht0=thenot*dr alt=alt_sfc call refracray(alt,tht0,r,rfndx,alt_101,n_101, &thtsfc,iapp,nplpts,layindx,losdist,rpl,mupl,x,ishadow,debug) amuo=cos(thtsfc) fs = amuo*z(ncp1) c write (6,'(i4,e12.4)') c &(i,-log(z(i))/exp(ttl(i)),i=1,ncp1) c c compute the energy input c engyin=1./omerf(sqchp*amuo,amuo,chpn) return end