subroutine firstz(h,ps,z,thenot,fracin,lmax,ncp1) c c*********************************************************************** c cccc c subroutine firstz 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. c 2.then using spline interpolation calculate the slant optical c paths in the model atmosphere. c c calling sequence- c call firstz(h,ps,z,thenot,fracin,fs,f1,f2,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 f1 r*8 o 3/16(cos(thenot)**2) c f2 r*8 o 3*sqrt(2)/16(sin(thenot)**2) 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) real *8 h(487),ps(487),zs(487),cofx(4,487) real*8 z(202) c include "consts.inc" include "out.inc" include "thkns.inc" include "chpmn.inc" 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) f1 = 0.1875*(1. + amuosq) f2 = 0.1875*sq2*(1. - amuosq) c c*****calculate slant opt. path of solar rays to reach c*****layer(i) of standard atmosphere c do 200 i = 1, lmax zs(i) = 0. raycon = sn*h(i) if (i .eq. 1) go to 190 duma= dsqrt((h(1) + raycon)*(h(1) - raycon)) c c*****add slant optical paths thru each layer(j) c*****lying between layer(1) and layer(i) c do 180 j = 2, i dumb= dsqrt((h(j) + raycon)*(h(j) - raycon)) path = fracin*(duma - dumb) zs(i) = zs(i) + path*dxs(j) 180 duma = dumb c c*****find the slant opt. path to reach layer 1 from the top of c*****the atmosphere- using chapman function approximation c 190 amu = raycon/h(1) amu = dsqrt(1. - amu**2) chpp = omerf(sqchp*amu, amu, chpn) chxx = omerf(sqchx*amu, amu, chxn) c c*****add the previous results to obtain total slant path zs to reach c*****each layer i. c zs(i) = zs(i) + chpp*ps(1) + chxx*(dxs(1) - ps(1)) 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 fs = amuo*z(ncp1) c write (6,'(i4,e12.4)') (i,z(i),i=1,ncp1) return end