c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/slant.f,v 2.22 2000/08/30 16:38:09 flittner Exp $ subroutine slant(h,ps,cofx,fracin,lmax) c c*********************************************************************** cccc c subroutine slant c c purpose- c 1. subdivide the layers of the standard atmosphere and c calculate the scattering and total optical depths of the c atmosphere from top to each layer. c 2. using the above results calculate the chapman constants. c c method- c 1. the layers are subdivided using exponential interpolation c 2. the total o.d. are calculated using spline interpolation c c calling sequence- c call slant(h,ps,cofx,fracin,juzsfr,lmax,layer) c c variable type i/o description c -------- ---- --- ----------- c c h(487) r*8 o height of the layers above surface(km) c height from earth center(fraction of r) c ps(487) r*8 o pressure at each layer (in atmosphere) c rayleigh opt. thickness of each layer c of the standard atmosphere c cofx(4,487) r*8 i spline interpolation coeff. c fracin o layer*radius of earth(km) c juzsfr i*4 i print flag c lmax i*4 i no. of layers in standard atmosphere c the following are input through contrl.cmn c layer i*4 i # of subdivisions to be made of each layer c of the standard atmosphere c hhold(101) r*8 i height of the layers above surface(km) height c from earth center(fraction of r) c pshold(101) r*8 i pressure at each layer (in atmosphere) rayleigh c opt. thickness of each layer of the standard c atmosphere c c other variables are returned thru common block 'chpmn' 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. Perform gravity correction when c computing rayleigh optical depth and storing in array ps. c Also use logical switch lgcorrect c to impliment the gravity correction to the rayleigh scattering c optical depth. c c last mod. Sep. 6, 1998...def c purpose: combine several commons into switches.cmn. c c last mod. Oct. 19, 1998...def c purpose: check to neg. and zero absorption optical depths in the c layers while computing xs and dxs (total extiniction values, which c should be at-least equal to the scattering optical depth.). c c last mod. May 6, 1999...def c purpose: pass pshold and hhold through contrl.cmn instead of commandline. c c last mod. Nov. 12, 1999...def c purpose: check for increase in absorber concentration at top of user profile c in calc. of the scale height. If it does, use scale height of cum. abs. c cccc implicit none c passed real*8 h(487),ps(487),cofx(4,487),fracin integer*4 lmax c commons include 'parameter.inc' include 'contrl.cmn' include 'atmos.cmn' include 'consts.cmn' include 'thkns.cmn' include 'chpmn.cmn' include 'switches.cmn' include 'atm_101.cmn' c real*8 interp_gc external interp_gc c local real*8 dum,hold,frac,dum1,dum2,dum3,holdc,holdd,scalp,scalx real*8 chp,chx,dabsmm,dabs2,scalxs integer*4 i,lmaxd2,j,k,mm,lmaxm1,layrm1,lmaxml c if(jprint(7).ne.0)write(33,6500) c c*****use the ps and h that were saved before c 80 do 82 i=1,lmax h(i)=hhold(i) 82 ps(i)=pshold(i) c c*****renumber the layers to put layer 1 near the top of the atmos. c***** (standard atmos. is read in a reverse order) c 85 continue lmaxd2 = lmax/2 do 100 i = 1, lmax if(gc_type.gt.0)then if(gc_type.eq.1)then ps(i) = beta*ps(i)*(1.0d0-pscaleforg*dlog(ps(i)))**2 else ps(i) = beta*ps(i)*interp_gc(ps(i)) endif else ps(i) = beta*ps(i) endif 100 continue c make sure that the value of dlog(ps(h=0))=tsl(101) ps(1)=exp(tsl(101)) if(jprint(7).ne.0)then write(33,*)'slant. ps after gc ' write(33,6450)(h(i),ps(i),i=1,lmax) write(33,*)dlog(ps(1)),tsl(101),dlog(ps(1))-tsl(101) endif do 110 i = 1, lmaxd2 k = lmax - i + 1 hold = h(i) h(i) = h(k) h(k) = hold hold = ps(i) ps(i) = ps(k) 110 ps(k) = hold c c*****subdivide the layers of std. atmos. by interpolation c*****assume pressure dependence between the layers of the form c***** p=pnot**(-h/hnot) c lmax = layer*(lmax - 1) + 1 lmaxm1 = lmax - 1 layrm1 = layer - 1 do 112 i = 1, lmax, layer j = lmax - i + 1 k = j/layer + 1 h(j) = h(k) 112 ps(j) = ps(k) frac = 1./float(layer) fracin = r*float(layer) lmaxml = lmax - layer do 114 i = 1, lmaxml, layer k = i + layer dum = (ps(k)/ps(i))**frac k = i do 114 j = 1, layrm1 k = k + 1 h(k) = h(k-1) - frac 114 ps(k) = ps(i)*dum**j c if(jprint(7).ne.0)then write(33,*)'slant. lmax now equals ',lmax write(33,'(1P2e12.4)')(h(i),ps(i),i=1,lmax) write(33,*)'slant. dlog(ps(lmax)) tsl(101)' write(33,*)dlog(ps(lmax)),tsl(101),dlog(ps(lmax))-tsl(101) endif c also make sure that xs(lmax) has a value initall equal to c tt(101) if(dlog(ps(lmax)).gt.tsl(101))then i = 101 k = i - 1 dum = tsl(i) dum1 = tsl(i) - dum dum2 = dum - tsl(k) dum3 = dum1*(cofx(1,k)*dum1**2 + cofx(3,k)) + dum2*(cofx(2,k)* 1dum2**2 + cofx(4,k)) xs(lmax) = dexp(dum3) endif c c*****use spline interpol. to obtain total opt. depth of each layer c***** in the standard atmosphere (xs) c j = 1 dum = dlog(ps(j)) do 130 i = 2, 101 if (dum .gt. tsl(i)) go to 130 k = i - 1 120 dum1 = tsl(i) - dum dum2 = dum - tsl(k) dum3 = dum1*(cofx(1,k)*dum1**2 + cofx(3,k)) + dum2*(cofx(2,k)* 1dum2**2 + cofx(4,k)) c check to see that xs is >= ps if(dum3.lt.dum)then xs(j)=dexp(dum) else xs(j) = dexp(dum3) endif j = j + 1 if (j .gt. lmax) go to 140 dum = dlog(ps(j)) if (dum .le. tsl(i)) go to 120 130 continue if(jprint(7).ne.0)then write(33,*)'slant. i j lmax xs',i,j,lmax,xs(lmax),xs(lmax-1) endif c c*****calculate ps(rayleigh opt. thickness) and dxs(total opt. thick) 140 dxs(1) = xs(1) holdc = ps(1) c calc. the scale height of the cum. abs. before ps is lost mm=20*layer+2 dabsmm=(xs(mm)-ps(mm)) dabs2=(xs(2)-ps(2)) if((dabs2.gt.0.0).and.(dabsmm.gt.0.0).and.(dabs2.lt.dabsmm))then scalxs=20./dlog(dabsmm/dabs2) else scalxs=0.0d0 endif c do 150 i = 2, lmax holdd = ps(i) ps(i) = holdd - holdc holdc = holdd dxs(i) = xs(i) - xs(i-1) c check to see that dxs(i) is ge ps(i) if(dxs(i).lt.ps(i))dxs(i)=ps(i) 150 continue do 160 i = 1, lmax xs(i) = dlog(xs(i)) 160 h(i) = 1. + h(i)*rinv if(jprint(7).ne.0) 1 write(33,6400)((h(i)-1.)*r,ps(i),xs(i),dxs(i),i=1,lmax) dum = r*h(1) c scale height for scattering c (use the change in the scattering optical c depth between layers mm and 2. We can't use layer 1 since it goes from c the pressure at the top of layer 2 to 0 at infinity scalp=20./dlog(ps(mm)/ps(2)) chp = dum/scalp*0.5 chpn = dsqrt(chp*pi) sqchp = dsqrt(chp) c scale height for absorption scalx = scalp dabsmm=(dxs(mm)-ps(mm)) dabs2=(dxs(2)-ps(2)) if((dabs2.gt.0.0).and.(dabsmm.gt.0.0).and.(dabs2.lt.dabsmm))then c the abs. optical depth is >0 in the layers and does decrease with alt. scalx=20./dlog(dabsmm/dabs2) chx = dum/scalx*0.5 chxn = dsqrt(chx*pi) sqchx = dsqrt(chx) if(jprint(7).ne.0)write(33,*)'slant. scale height o3 ', &'between layer mm and 2 is ',scalx else c must do somethingelse c simplest is set to zero scalx=0.0d0 chx=0.0d0 chxn=0.0d0 sqchx=0.0d0 dum2=dxs(1)-ps(1) if(dum2.gt.0.0)then c however, if dxs(1)-ps(1) is greater than zero, we should calc. something c So, will average over the last 20 km. j=0 dum2=dxs(2)-ps(2) do i=3,mm dum3=dxs(i)-ps(i) if((dum2.gt.0.0).and.(dum3.gt.0.0).and.(dum3.gt.dum2))then j=j+1 scalx=r*(h(i-1)-h(i))/log(dum3/dum2)+scalx endif dum2=dum3 enddo if(j.gt.0)then scalx=scalx/dble(j) chx = dum/scalx*0.5 chxn = dsqrt(chx*pi) sqchx = dsqrt(chx) if(jprint(7).ne.0)write(33,*)'slant. scale height o3', &' averaged over the last 20 km is ',scalx else c will use the scale height for the cum. abs. if(scalxs.gt.0.0)then chx = dum/scalxs*0.5 chxn = dsqrt(chx*pi) sqchx = dsqrt(chx) if(jprint(7).ne.0)write(33,*)'slant. scale height', &' o3 for cum. abs. is ',scalxs endif endif endif endif if(jprint(7).ne.0) write(33,*)'slant. scalp scalx',scalp,scalx 6400 format (1h1,2(5x,1hh,14x,2hps,14x,2hxs,13x,3hdxs,9x)/1h,/, 1 (1h ,2(1Pe10.4, 3d16.4, 5x))) 6450 format (1P2e14.4) 6500 format('slant. debugging output-------') 6550 format(i9,3E12.4) return end