ciniSegAtoms.f

Go to the documentation of this file.

      subroutine cinisegatmos
!        read the atmosphere data
!
      implicit none
#include  "Zmanager.h"
#include  "Zmanagerp.h"
#include  "Zearth.h"
#include  "Zatmos.h"

      integer icon,  ios, i
      character*100 msg

      call copenf(tempdev, atmosfile, icon)
      if(icon .ne. 0) stop 9999
      call cskipcomment(tempdev, icon)
      if(icon .ne. 0) stop 9999

      nodes = 0
      mostz = 0
      do while( .true. )
         read(tempdev, *, iostat=ios)
     *      atmos%z(nodes+1), atmos%T(nodes+1), atmos%P(nodes+1), 
     *      atmos%rho(nodes+1),  atmos%a(nodes+1), atmos%d0(nodes+1),
     *      atmos%cumd(nodes+1),  atmos%H(nodes+1)
         if(ios .ne. 0) goto 10
         if(nodes .ge. maxnodes) then
            write(msg,*) 'numbr of nodes for atmosphere > ', maxnodes
            call cerrormsg(msg, 0)
         endif
         nodes = nodes + 1
         if(mostz .eq. 0 .and.   znode(nodes) .gt. 2900.) then
            mostz = nodes
         endif
!        write(0,*) Znode(nodes), Anode(nodes)
      enddo
 10   continue
!     write(*, *) " # of nodal points =", nodes, ' mostz', mostz

      do i=1, nodes-1
         if(anode(i) .ne. 0.) then
            fd1i(i) = fd1(znode(i+1), anode(i), znode(i), 
     *           scalehnode(i))
            rhopnode(i) =  rhonode(i) *(-1.d0 -1.d0/anode(i)) *
     *           anode(i)/scalehnode(i) 
            pwp(i) =-2.d0-1.d0/anode(i)
         else
            fd0i(i) = fd0( znode(i+1), znode(i), scalehnode(i) )
         endif
      enddo
      endif
      end
!     *********************************   
      real*8 function  cvh2den(z)
!     *********************************   
!         gives density of air at a given vertical height.

!
      implicit none
#include  "Zatmos.h"

      real*8 z     ! input.  vertical height in m
!                   function value:  density in kg/m3.
      integer i 

      if( z .ge. znode(mostz) .and. z .lt. znode(mostz+1)) then
         i = mostz
      elseif( z .ge znode(mostz + 1) .and. z .lt. znode(mostz+2) then
         i = mostz + 1
      else
         do i = mostz, 2, -1
            if(z .ge. znode(i-1) .and. z .lt. znode(i)) then
               goto 100
            endif
         enddo
         call kdwhereis(z, nodes, znode, 1, i)
         if(i .eq. 0) then
            call cerrormsg('height becomes too small ', 0)
         endif
      endif
 100  continue
      if(anode(i) .ne. 0.)then
         cvh2den = 
     *        rhonode(i) * 
     *       (1.0 + anode(i)* 
     *          (z-znode(i))/scalehnode(i) )**(-1.0-1.d0/anode(i)) 
      else
         cvh2den =
     *        rhonode(i) * exp(- (z-znode(i))/scalehnode(i))         
      endif
      end
!     *************************
      real*8 function cvh2thick(z)
!           vertical height to thickness of air above that point.
!      z:  input. vertical height in m
!  function value. output.  thickness in kg/m2.
!
!
      implicit none
#include  "Zearth.h"
#include  "Zstdatmos.h"
#include  "Zatmos.h"

      real*8 z        ! input
      integer i

#include  "Zstdatmosf.h"


!          The gramage between  given heights, z1 and z2  is by
!
!              d = D0node *(fd(z1) - fd(z2))  where
!         
!            fd(z) = (1+ a(z-z0)/H(z0))**(-1/a)                   (a != 0)
!
!                  =  exp(-(z-z0)/H )                             (a = 0)

      if( z .ge. znode(mostz) .and. z .lt. znode(mostz+1)) then
         i = mostz
      elseif( z .ge znode(mostz + 1) .and. z .lt. znode(mostz+2) then
         i = mostz + 1
      else
         do i = mostz, 2, -1
            if(z .ge. znode(i-1) .and. z .lt. znode(i)) then
               goto 100
            endif
         enddo
         call kdwhereis(z, nodes, znode, 1, i)
         if(i .eq. 0) then
            call cerrormsg('height becomes too small ', 0)
         endif
      endif
 100  continue

      if(anode(i) .ne. 0.)then
         cvh2thick = dsumnode(i) -
     *     d0node(i) * (
     *     1.0- fd1(z, anode(i), znode(i), scalehnode(i)) 
     *                 )
      else
         cvh2thick = dsumnode(i) +
     *     d0node(i) * (1. -
     *       fd0(z, znode(i), scalehnode(i) )
     *                 )
      endif
      end
!     *************************
      real*8 function cvthick2h(t)
!         vertical thickness to vertical height conversion
!     t:  input. vertical tikness in kg/m2
!  function value. output.  height in m.
!
!
      implicit none
#include  "Zatmos.h"

      integer i

      real*8  t,  temp


      if( z .ge. znode(mostz) .and. z .lt. znode(mostz+1)) then
         i = mostz
      elseif( z .ge znode(mostz + 1) .and. z .lt. znode(mostz+2) then
         i = mostz + 1
      else
         do i = mostz, 2, -1
            if(z .ge. znode(i-1) .and. z .lt. znode(i)) then
               goto 100
            endif
         enddo
         call kdwhereis(z, nodes, znode, 1, i)
         if(i .eq. 0) then
            call cerrormsg('height becomes too small ', 0)
         endif
      endif
 100  continue

      if(anode(i) .ne. 0. )then
!         cvh2thick = Dsumnode(i) -
!     *     D0node(i) * (
!     *     1.0- fd1(z, Anode(i), Znode(i), Scalehnode(i)) 
!     *                 )
!                    solve above eq.
         temp =1.0 -  (dsumnode(i) - t )/d0node(i)

         cvthick2h = (1.0 -  temp**(-anode(i)))*
     *               scalehnode(i)/anode(i) + znode(i)
      else
!          t = Dsumnode(i) +
!     *     D0node(i) * (1. -
!     *       fD0(z, Znode(i), Scalehnode(i) )
!     *                 )
!            solve above eq.
         temp =1.0- (t -dsumnode(i)) /d0node(i) 
!                   temp =  exp(-(z-z0)/H )    
         cvthick2h =znode(i)- log(temp)*scalehnode(i)
      endif
      end
!     *********************************   
      real*8 function  cvh2denp(z)
!     *********************************   
!         gives derivative of the density of
!     air at a given vertical height.
!      h: vertical height in m
!      function value: dRhonode/dz density in kg/m4
!
      implicit none
c----      include 'Zearth.h'
#include  "Zearth.h"
c----      include 'Zstdatmos.h'
#include  "Zstdatmos.h"
!
      integer i

      real*8 z


      if( z .ge. znode(mostz) .and. z .lt. znode(mostz+1)) then
         i = mostz
      elseif( z .ge znode(mostz + 1) .and. z .lt. znode(mostz+2) then
         i = mostz + 1
      else
         do i = mostz, 2, -1
            if(z .ge. znode(i-1) .and. z .lt. znode(i)) then
               goto 100
            endif
         enddo
         call kdwhereis(z, nodes, znode, 1, i)
         if(i .eq. 0) then
            call cerrormsg('height becomes too small ', 0)
         endif
      endif
 100  continue

      if(anode(i) .ne. 0.)then
         cvh2denp = 
!       *             Rhonode(i-1) *(-1.D0node -1.D0node/Anode(i-1)) *
!       *             Anode(i-1)/Scalehnode(i-1) *     = Rhonodep(i)
     *        rhopnode(i) *
     *        (1.d0 + anode(i)* 
!     *             (z-Znode(i))/Scalehnode(i) )**(-2.d0-1.d0/Anode(i)) 
     *       (z-znode(i))/scalehnode(i) )**pwp(i)
      else
         cvh2denp =
     *        -rhonode(i) * exp(- (z-znode(i))/scalehnode(i))         
     *        /scalehnode(i)
      endif
      end
!     *********************************   
      real*8 function  cvh2den2p(z)
!     *********************************   
!         gives double derivative of the density of
!     air at a given vertical height.
!      h: vertical height in m
!      function value: d (dRhonode/dz)/dz density in kg/m5
!
      implicit none
c----      include 'Zearth.h'
#include  "Zearth.h"
c----      include 'Zstdatmos.h'
#include  "Zstdatmos.h"
!
      integer i

      real*8 z

      if(first) call cstdatmos0

!          check if z < 11.2km
      if(z .lt. znode(2)) then
         cvh2den2p = rhonode00 * pw*(pw-1.d0node)/hl/hl
     *          * ( (ha - z)/hl )**(pw-2.d0node)
      else
         do i = 3, nodes
            if(z .lt. znode(i) .or. i .eq. nodes ) then
               if(anode(i-1) .ne. 0.)then
                  cvh2den2p = 
     *               rhonodep(i) * pwp(i) *anode(i-1)/scalehnode(i-1)*
     *             (1.d0node + anode(i-1)* 
     *             (z-znode(i-1))/scalehnode(i-1) )**(pwp(i)-1.d0node)
                   goto 10
               else
                  cvh2den2p =
     *            rhonode(i-1) * exp(- (z-znode(i-1))/scalehnode(i-1))         
     *            /scalehnode(i-1)**2
                  goto 10
               endif
            endif
         enddo
 10      continue
      endif
      end