catmosCnst1.f

Go to the documentation of this file.

      subroutine catmoscnst1

!   catmosCnst1:  compute some basic constants
!      a,  b,  d0, cum, H
! where
!     a, b; const shown below
!       d0: see the formula below.
!     cumd: amount of atmospher above the node.
!        H: scale height at the node
!     all in mks unit.
!
!     The scale height is approximated by a
!   number of stright lines as a function of height. The data in
!   stdatmos1.d gives height, temperatur, etc at  each nodal point.
!
!         The scale height, H, is expressed by H = H0 + a(z-z0)
!                                                = kT/mg      
!   in each region.
!   We neglect height dependence of gravitational accelleration g,
!   and the average mass of  air molecules, m.   
!   Since the data table gives T(z)= T0 + b(z-z0) at the nodal points, 
!   we can first get b,
!   and then a by a = dH/dz = k/mg * b.  H at a nodal point, z,  
!   is obtained as H(z) =kT(z)/mg.
!
!   The density off a nodal point is given by
!
!              rho = rho0 * (1+ a(z-z0)/H(z0))**(-1-1/a)      (a != 0)
!                  = rho0 * exp(- (z-z0)/H)            (a =0; hence H is const)
!
!   (We employ H(z0) as the scale height in the segment)
!   The amount of air between  given heights, z1 and z2  is by
!
!              d = d0 *(fd(z1) - fd(z2))  where
!         
!            fd(z) = (1+ a(z-z0)/H(z0))**(-1/a)                   (a != 0)
!
!                  =  exp(-(z-z0)/H )                             (a = 0)
!
!   where  d0 = rho0*H(z0)
!
!   If z1=z0,  d becomes
!
!              d= d0 ( 1 - fd(z2))
!
!

!
!    
      implicit none

#include  "Zglobalc.h"
#include  "Zatmos.h"

      integer i

      

      real*8 m, g, k, nuc, c2, kbymg, dtdz
!      m: avrage molecule mass number of air
!      g: gravitational acceleration.
!      k: Boltzman's const.
!     c2: c**2
      parameter( m=14.5d0 *2, nuc =938.3d6, g=9.80665d0, c2=c*c,
     *          k= 8.617d-5, kbymg = k*c2/(m*nuc*g))

      integer nodes

#include  "Zstdatmosf.h"

      nodes = atmos%nodes
      do i = 1, nodes-1
         dtdz = (atmos%T(i+1) - atmos%T(i))/ (atmos%z(i+1) - atmos%z(i))
         atmos%H(i) = kbymg * atmos%T(i)                 ! at atmos%z(i)
         atmos%b(i) = dtdz
         atmos%a(i) = kbymg * atmos%b(i)
         atmos%d0(i) = atmos%rho(i) * atmos%H(i)
!///////////
!        write(0, *)
!     * atmos.z(i), atmos.d0(i), atmos.a(i), atmos.b(i), atmos.H(i)
!////////
      enddo

      atmos%cumd(nodes) =  atmos%rho(nodes) * hinf   ! put  very small amout

      do i = nodes-1, 1, -1
         if(atmos%a(i) .eq. 0.) then
!              exponential atmosphere
            atmos%cumd(i) = atmos%cumd(i+1) +
     *       atmos%d0(i)*
!     *       (1.0- exp(-(atmos.z(i+1) - atmos.z(i))/atmos.H(i)))
     *       (1.0- fd0(atmos%z(i+1), atmos%z(i), atmos%H(i)))
         else
            atmos%cumd(i) = atmos%cumd(i+1) +
     *        atmos%d0(i)*
     *        (1.0 -
     *          fd1(atmos%z(i+1), atmos%a(i),atmos%z(i), atmos%H(i) ) )

!     *        (1.0+atmos.a(i)*(atmos.z(i+1) -atmos.z(i))/atmos.H(i))
!     *        **(-1./atmos.a(i))
!     *         )
         endif
      enddo
      do i = 1, nodes
!         atmos.logP(i) = log(atmos.P(i))
         atmos%logrho(i) = log(atmos%rho(i))
         atmos%logcumd(i) = log(atmos%cumd(i))
         atmos%logcumdi(atmos%nodes-i+1) = atmos%logcumd(i)
         atmos%logrhoi(atmos%nodes-i+1) = atmos%logrho(i)
         atmos%zi(atmos%nodes-i+1) = atmos%z(i)
      enddo
      end

      subroutine catmoscnst2
!         compute c-spline coef. for later use
      implicit none
#include "Zatmos.h"

      integer  nodes
      
      nodes = atmos%nodes
!                           height--> rho
      call kcsplcoef(atmos%z, atmos%logrho, nodes, atmos%coefh2r,
     *      maxnodes)
!                           height--> log(depth)
      call kcsplcoef(atmos%z, atmos%logcumd, nodes, atmos%coefh2d,
     *      maxnodes)
!                           h--> log(P)
!      call kcsplCoef(atmos.z, atmos.logP, nodes, atmos.coefh2P,
!     *      maxnodes)
!                           h--> scale H
      call kcsplcoef(atmos%z, atmos%H, nodes-1, atmos%coefh2H,
     *      maxnodes)
!                           h--> T
      call kcsplcoef(atmos%z, atmos%T, nodes, atmos%coefh2T,
     *      maxnodes)
!                           log(rho)--> h
!      call kcsplCoef(atmos.logrho, atmos.z, nodes, atmos.coefr2h,
!     *      maxnodes)
!                           log(depth)--> h
      call kcsplcoef(atmos%logcumdi, atmos%zi, nodes, atmos%coefd2h,
     *      maxnodes)

!               log(detph) --> log(rho)
      call  kcsplcoef(atmos%logcumdi, atmos%logrhoi, nodes,
     *      atmos%coefd2r, maxnodes)
!                           log(P) --> h 
!      call kcsplCoef(atmos.logP, atmos.z, nodes, atmos.coefP2h,
!     *      maxnodes)

      end