c*****************************************************
c* extracts O and N2 mass denisties a *.egw file     *
c* ( EB / January 2026 )                             *
c*****************************************************
c
      implicit none
      integer lev1,nlon,nlat,nho,idat
      parameter(lev1=261,nlon=360,nlat=180,nho=nlon*nlat,
     &          idat=16  )
      integer i,j,l,it,id,krec,l1,l2
      real*4  dat(nlon,nlat,lev1,idat),alt(lev1)
      character ch5*5(37) 
c      
      data ch5 /'12:00','12:10','12:20','12:30','12:40','12:50',
     &          '13:00','13:10','13:20','13:30','13:40','13:50',
     &          '14:00','14:10','14:20','14:30','14:40','14:50',
     &          '15:00','15:10','15:20','15:30','15:40','15:50',
     &          '16:00','16:10','16:20','16:30','16:40','16:50',
     &          '17:00','17:10','17:20','17:30','17:40','17:50',
     &          '18:00'/
      do i=1,37
        print*,'i,ch5=',i,ch5(i)
      enddo
c
c altitudes (km) 
c     delta z = 0.5 km up to 70 km
      do l = 1,140
        alt(l) = 0.5d0*dble( (l-1) )
      enddo
c    delta z = 1 km from 70 to 120 km
      do l = 1,50
        alt(140+l) = 70d0  + dble( (l-1) )
      enddo
c    delta z = 2 km from 120 to 200 km
      do l = 1,40
        alt(190+l) = 120d0 + dble( 2*(l-1) )
      enddo
c    delta z = 5 km from 200 to 300 km
      do l = 1,20
        alt(230+l) = 200d0 + dble( 5*(l-1) )
      enddo
c     delta z = 10 km from 300 to 400 km
      do l = 1,11
        alt(250+l) = 300d0 + dble( 10*(l-1) )
      enddo
      print*,'alt = ',real(alt)
c
      l1 = 196 ! 130 km
      l2 = 221 ! 180 km      
      print*,'alt(l1) = ',alt(l1) 
      print*,'alt(l2) = ',alt(l2) 
c      
      open(12,file='29DEC2020.12-18.20.10m.sz.grads',
     &        form='unformatted',access='direct',
     &        recl=nlon*nlat, status='old') 
c
      krec = 0
      DO IT = 13,19 ! from 12 UT to 15 UT
c
      open(10,file='29DEC2020.12-18.chem.'//ch5(it)//'UT.txt',
     &        status='unknown') 
      write(10,'(1x,i2,a)')l2-l1-1,' altitudes (km) ='
      write(10,*)( int(alt(l)), l=l1,l2 ) 
      write(10,*)'180 latitudes (deg) = '   
      write(10,*)( -89.5 + real(j-1) , j=1,nlat ) 
      write(10,*)'360 longitudes (deg) = '   
      write(10,*)(i, i=1,nlon) 
       
        do id = 1,idat
        do  l = 1,lev1 
          krec = krec+1
          read (12,rec=krec)
     &      ( ( dat(i,j,l,id),i=1,nlon ) , j=1,nlat )
        enddo; enddo
        print*,'after reading time stamp, it = ',it
c
c       make printout of mass densities 
c       in 10**(-9) kg/m**3 = 10**(-12) g/cm**3 
c       from 130 to 180 km (l1 to l2)
c       (1) O2
        write(10,*)'O2 mass density' 
        do l = l1,l2
          write(10,*) 
     &    ( ( 1e9*dat(i,j,l,6)*dat(i,j,l,7), i=1,nlon ), j=1,nlat ) 
        enddo 
c       (2) O
        write(10,*)'O mass density' 
        do l = l1,l2
          write(10,*) 
     &    ( ( 1e9*dat(i,j,l,6)*dat(i,j,l,8), i=1,nlon ), j=1,nlat ) 
        enddo 
c       (3) N2 
        write(10,*)'N2 mass density'
        do l = l1,l2
          write(10,*)
     &    ( ( 1e9*dat(i,j,l,6)*(1.0-dat(i,j,l,7)-dat(i,j,l,8)), 
     &    i=1,nlon ), j=1,nlat )
        enddo
        close(10) 
      ENDDO 
      stop'ok'
      end
