spathz.f

Go to the documentation of this file.

      subroutine spathz(amu0,htp,taur,taum,taua,totsp,nmodl)            
c                                                                       
      implicit real*8 (a-h,o-z)                                         
      real*8 htp(1000),taur(1000),taum(1000),taua(1000),                
     1       totsp(1000),fxsp(1000)                                     
      real*8 dht(1000),dtr(1000),dtm(1000),dtot(1000),sum(1000)         
      real*8 hs(1000),dta(1000),yy,xx
c                                                                       
      pi=dacos(-1.0d0)                                                  
      conv=pi/180.0d0                                                   
      snmu0=dsqrt(1.0-amu0**2)                                          
c     scale heights                                                     
      shp=7.930024d0                                                    
      sho=4.178129d0                                                    
      re=6371.0d0                                                       
c     compute (h+re)/re for each pressure level and geometric thickness 
c     of each layer                                                     
      do i=1,nmodl                                                      
         hs(i)=(re+htp(i))/re                                           
         if(i.lt.nmodl)then                                             
            dht(i)=htp(i)-htp(i+1)                                      
            dtr(i)=taur(i+1)-taur(i)                                    
            dtm(i)=taum(i+1)-taum(i)                                    
            dta(i)=taua(i+1)-taua(i)                                    
         endif                                                          
c     write(6,801)i,dht(i),dtr(i),dtm(i),dta(i)                         
801   format('i,dht,dtr,dtm,dta',i5,1p4e12.4)                           
      enddo                                                             
c     compute the slant path                                            
c                                                                       
      sdp=(re+htp(1))/shp                                               
      sdo=(re+htp(1))/sho                                               
      do i=1,nmodl                                                      
         sum(i)=0.0d0                                                   
         aa=snmu0*hs(i)                                                 
           if(i.eq.1)then                                               
             call chapmn(sdp,amu0,chpp)                                 
             call chapmn(sdo,amu0,chpo)                                 
             sum(i)=sum(i)+taur(1)*chpp+taua(1)*chpo                    
             totsp(i)=sum(i)                                            
             fxsp(i)=amu0*exp(-sum(i))                                  
          else                                                          
             xx=sqrt((hs(1)+aa)*(hs(1)-aa))                             
             do j=2,i                                                   
                yy=sqrt((hs(j)+aa)*(hs(j)-aa))                          
                ds=(xx-yy)*re                                           
                dst=ds/dht(j-1)                                         
                sum(i)=sum(i)+dst*(dtr(j-1)+dtm(j-1)+dta(j-1))          
c      write(6,803)i,j,xx,yy,ds,dst,sum(i)                              
803    format('i,j,xx,yy,ds,dst,sum',2i4,1p5e10.2)                      
                xx=yy                                                   
             enddo                                                      
               snmu=aa/hs(1)                                            
               amu=dsqrt(1.0-snmu**2)                                   
               call chapmn(sdp,amu,chpp)                                
               call chapmn(sdo,amu,chpo)                                
               sum(i)=sum(i)+taur(1)*chpp+taua(1)*chpo                  
               totsp(i)=sum(i)                                          
               fxsp(i)=amu0*dexp(-sum(i))                               
c      write(6,802)i,snmu,amu,chpp,sum(i),fxsp(i)                       
802    format('i,snmu,amu,chpp,sum(i),fxsp(i)',                         
     1          i5,1p5e11.3)                                            
c      write(6,809)i,sum(i),htp(i),fxsp(i)                              
809    format('i,sum,htp,fxsp',i3,1x,1p3e12.4)                          
          endif                                                         
      enddo                                                             
      return                                                            
      end                                                               
c***********************************************************************