      program rdthin

!     Add ice and land    Michiko Masutani 7/29/2009
!     This program read NR Thin data. The code combineis rdthin.f and
!     rdprof60agl2.f90.  - Tong 3/10/2008

!     This program read NR DBL91 data, with GOSE obs's lon, lat, hour
!     - Start Tong Zhu 11/08/2007
!     Add compute saza, soza - Tong Zhu 11/16/2007
!     Add time,cloud_Total,cloudL,cloudM,cloudH - Tong Zhu 12/04/2007

      implicit none
      integer n, l, k, m, i
      integer nclm, nsfc, nupr, kmax, nsat, lugrb

      integer msfc, msat, isvar
      parameter (msfc=100,msat=50)

      character(80) filenm,istxt
      common /sfcobs/ nsfc,isvar(msfc),istxt(msfc)
      character(50) satnd(16),climd(4),surfd(11)
      character(10) upprd(6),upprv(6)
!      dimension clim(4),satn(msat),surf(msfc),psig(91),uppr(91,6)
      real*4 clim(4),satn(msat),surf(msfc),psig(91),uppr(91,6)
      data lugrb/10/

      real gg
      data gg /9.80665/     !m/s**2 gravity constant
      real plev(91), ptop, ttop
      real terr,lat,lon,hour,u10m,v10m
      real ismk,almk,sktp,snow(6)
! 
!   Variables used by GMAO 
!......................................
! ID /31,142,143,165,166,167,168,172,186,187,188,235/
!      1,17, 18, 28, 29, 30, 31, 33, 42, 43, 44, 57
! var !haracters for different fields:
!  ismk: surface sea-ice cover (0,1)
!  rain: surface stratiform precip. (m)
!  conp: surface convective precip. (m)
!  u10m: uwind at 10 meter
!  v10m: vwind at 10 meter
!  tp2m: surface temperature at 2 meter
!  dt2m: dew point temperature at 2 meter
!  almk:  sea (Aqua) land mask (0,1)
!  lcld: low cloud cover fraction (0,1)
!  mcld: middle  cloud cover fraction (0,1)
!  hcld: high cloud cover fraction (0,1)
!  sktp:  surface skin temperature  (k)
!......................................
!  snow:
!  1   2   kpds=032    Snow albedo [(0-1)]
!  2   3   kpds=033    Snow density [kg m**-3]
!  3   5   kpds=044    Snow evaporation [m of water
!  4   6   kpds=045    Snowmelt [m of water]
!  5  16   kpds=141    Snow depth [m of water equiv
!  6  19   kpds=144    Snowfall (convective + strat
!
     
      real lveg_c,hveg_c,lveg_t,hveg_t
      real saza,soza,sazm,cldt,cldl,cldm,cldh

      integer year, month, day, jday   !Julian Day
      integer n_prof,n_aero,nct

      character(6) sate
      include 'satellite.inc'
      real satlon, satlat

!  define the satellite elements
!  -----------------------------

      satnd( 1) = '004001  YEAR  YEAR         YEAR                   '
      satnd( 2) = '004002  MNTH  MONTH        MONTH                  '
      satnd( 3) = '004003  DAYS  DAY          DAY                    '
      satnd( 4) = '004004  HOUR  HOUR         HOUR                   '
      satnd( 5) = '004005  MINU  MINUTE       MINUTE                 '
      satnd( 6) = '004006  SECO  SECOND       SECOND                 '
      satnd( 7) = '006002  CLON  DEGREES      LONGITUDE              '
      satnd( 8) = '005002  CLAT  DEGREES      LATITUDE               '
      satnd( 9) = '001007  SAID  CODE TABLE   SAT IDENTIFIER         '
      satnd(10) = '002019  SIID  CODE TABLE   SAT INSTRUMENTS        '
      satnd(11) = '005043  FOVN  NUMERIC      BEAM POSITION          '
      satnd(12) = '008012  LSQL  CODE TABLE   LAND/SEA QUALIFIER     '
      satnd(13) = '007024  SAZA  DEGREE       SAT ZENITH ANGLE       '
      satnd(14) = '007025  SOZA  DEGREE       SOLAR ZENITH ANGLE     '
      satnd(15) = '010001  HOLS  METER        HEIGHT OF LAND SURFACE '
      satnd(16) = '007002  HMSL  METER        HEIGHT OR ALTITUDE     '

!  define the climate values
!  -------------------------

      climd(1)='iv=27 ! low  vegetation cover'
      climd(2)='iv=28 ! high vegetation cover'
      climd(3)='iv=29 ! low  vegetation type'
      climd(4)='iv=30 ! high vegetation type'

!  read in the list of surface values to include
!  ---------------------------------------------

      open(lugrb,file='surface.lst');n=0
!!!3     if(n+1.gt.msfc) call bort('too many surface definitions')
3     if(n+1.gt.msfc) stop
      read(lugrb,'(5x,i3,4x,a80)',end=4)isvar(n+1),istxt(n+1)
      n=n+1; nsfc=n
      goto 3
4     print*,'read ',nsfc,' surface variable definitions'
!!!      if(nsfc.gt.msfc) call bort('too many surface definitions')
      if(nsfc.gt.msfc) stop

!  define the upper air values
!  ---------------------------

      upprd(1)=' cloud cov'; upprv(1)='    iv=248'
      upprd(2)=' cloud ice'; upprv(2)='    iv=247'
      upprd(3)=' cloud h2o'; upprv(3)='    iv=246'
      upprd(4)=' ozone mmr'; upprv(4)='    iv=203'
      upprd(5)=' tempature'; upprv(5)='    iv=130'
      upprd(6)=' spf humid'; upprv(6)='    iv=133'

!  open the thinned dbl91 file and read one record
!  -----------------------------------------------

!!!      read(5,'(a)') file; print*; print*,'reading ',file(1:60)
!!!      open(20,file=file,form='unformatted')
!     filenm='noaa15_amsua_2005050200.dbl91'
!     filenm='goes12_sndr_2005050200.dbl91'
!      filenm='goes.thin.dbl91'
!      open(20,file=filenm,form='unformatted')
      open(20,form='unformatted')
      open(12,file='nr2goes.thin.input',status='unknown',form='formatted')
      open(13,file='ptsN46L.txt',status='unknown',form='formatted')

      n_aero=0
      n_prof=0; nct=0

      do n=1,403239
!     do n=1,3

      read(20,end=999) nsat,nclm,nsfc,nupr,kmax,(satn(m),m=1,nsat), &
                    (clim(m),m=1,nclm),(surf(m),m=1,nsfc), &
                    (psig(l),l=1,kmax),((uppr(l,m),l=1,kmax),m=1,nupr)

      n_prof=n_prof+1

      year=satn(1)
      month=satn(2)
      day=satn(3)
      hour=satn(4)+satn(5)/60.+satn(6)/3600.

      lon=satn(7)
      lat=satn(8)

      saza=satn(13)
      soza=satn(14)
      sazm=0.0
!      added by msq
      lveg_c=clim(1)   
      hveg_c=clim(2)
      lveg_t=clim(3)
      hveg_t=clim(4)

      if(sate.eq.'goes10'.or.sate.eq.'goes12') then
         if(n_prof.eq.1) then
            call date2jday(year,month,day,jday)
            print*,'Process ',sate,' for date:',year,month,day,hour,' JDay:',jday
         endif

         if(sate.eq.'goes10') then
           satlon=-135.0
           satlat=0.0
         endif
         if(sate.eq.'goes12') then
           satlon=-75.0
           satlat=0.0
         endif

         call sazasoza(jday,satlon,satlat,lon,lat,hour,saza,soza)

      endif

      ttop=uppr(1,5)

      terr=surf(13)/gg
      u10m=surf(28)
      v10m=surf(29)
      cldt=surf(27)
      cldl=surf(42)
      cldm=surf(43)
      cldh=surf(44)
      ismk=surf(1)   ! msq 090729
      almk=surf(33) ! msq 090729
      sktp=surf(56)
      snow(1)=surf(2)
      snow(2)=surf(3)
      snow(4)=surf(5)
      snow(5)=surf(16)
      snow(6)=surf(19)


      do l=1,kmax
         psig(l)=psig(l)/100.0
      enddo
      do l=1,kmax-1
         plev(l)=(psig(l)+psig(l+1))/2.0
      enddo
      plev(kmax)=surf(25)/100.0

      ptop=0.005   !hPa       (Liu)
      
!!!     print*,'nsat,nclm,nsfc,nupr,kmax=',nsat,nclm,nsfc,nupr,kmax
!!!     print*
!!!     print'(f20.2,2x,a50)',(satn(m),satnd(m),m=1,16)                   
!!!     print*
!!!     print'(f20.5,2x,a50)',(clim(m),climd(m),m=1,nclm)
!!!     print*,'----------'
!!!     print'(f20.5,2x,2i5,2x,a50)',(surf(m),isvar(m),m,istxt(m),m=1,nsfc)
!!!     print*,'----------'
!!!9     format(f20.5,2x,2i5,2x,a50)
!!!     print*
!!!     print'("    pres(pa)",6a12)',upprd
!!!     print'("            ",6a12)',upprv
!!!     print*
!!!     print'(7e12.4)',(psig(k),(uppr(k,m),m=1,nupr),k=1,kmax,10)
!!!     print*
      
!      print*,n_prof,terr,lat,lon,psig(1),ttop,u10m,v10m
!      print*,saza,soza,sazm,hour,cldt,cldl,cldm,cldh
!      print*,n_prof,lon,lat,hour,saza,soza

      if(n_prof.eq.1) then
        print*,'n_prof=',n_prof,year,month,day,hour,lat,lon,saza,soza
      endif
      write(12,'(4i12)') year, month, day, n_aero
      write(12,'(15f12.4)') terr,lat,lon,psig(1),ttop,u10m,v10m,&
                            saza,soza,sazm,hour,cldt,cldl,cldm,cldh
      write(12,'(9f12.4)')ismk,almk,sktp,(snow(i),i=1,6)
      write(12,'(4f12.4)')lveg_c,hveg_c,lveg_t,hveg_t
      do k=1,kmax,2
      write(12,'(7e13.4)') plev(k),psig(k),uppr(k,5),uppr(k,6),uppr(k,4),&
                           uppr(k,3),uppr(k,2)
      enddo

      enddo        !------ end of n=1,max-profiles ------
999   print*,'n_prof=',n_prof,year,month,day,hour,lat,lon,saza,soza

      write(13,'(1x,i12,1x,i12)') 46, n_prof

      close(10)
      close(12)
      close(13)
      close(20)

10    stop
      end
!
