      program read radiation

      implicit none

      integer                  nx    ,ny
      parameter               (nx=250,ny=250)
      real      shortwave_down(nx    ,ny    ),
     &          surf_alb      (nx    ,ny    ),
     &          par_down      (nx    ,ny    ),
     &          surf_alb_v    (nx    ,ny    ),
     &          longwave_net  (nx    ,ny    ),
     &          rflvs         (nx    ,ny    ),
     &          toa_swd       (nx    ,ny    ),
     &          cloud_alb     (nx    ,ny    ),
     &          surf_temp     (nx    ,ny    ),
     &          pw            (nx    ,ny    ),
     &          rfclr         (nx    ,ny    ),
     &          lat           (nx    ,ny    ),
     &          lon           (nx    ,ny    ),

c      Derived variables.

     &          shortwave_up  (nx    ,ny    ),
     &          shortwave_net (nx    ,ny    ),
     &          par_up        (nx    ,ny    ),
     &          par_net       (nx    ,ny    ),
     &          longwave_up   (nx    ,ny    ),
     &          longwave_down (nx    ,ny    ),
     &          net_rad       (nx    ,ny    )

      real      scale         (11           )

      real      eps,sigma
      integer   i,j,julianday,hour,min
      character fname*60

c      Integer conversion scale factor for each variable.

      data scale/0.0500,0.0001,0.0500,0.0001,0.0500,
     &    0.0001,0.0500,0.0001,0.0100,0.0100,0.0001/

c      Read the lat/lon navigation data.

      open(1, file  ='latlon.dat',
     &        form  ='formatted' ,
     &        status='old'       )
      read(1,'(10(f7.3,f9.3))') ((lat(i,j),lon(i,j),i=1,nx),j=1,ny)
      close (1,status='keep')

c      Open the direct access file. 

      call getarg(1,fname)      
      open (2,file  =fname        ,
     &        status='old'        ,
     &        form  ='unformatted',
     &        access='direct'     ,
     &        recl  =2*nx*ny/10   )

c      Read the header to get the date and time.

      call header(julianday     ,hour     ,min)
      call read1 (shortwave_down,scale( 1),  2)
      call read1 (surf_alb      ,scale( 2), 12)
      call read1 (par_down      ,scale( 3), 22)
      call read1 (surf_alb_v    ,scale( 4), 32)
      call read1 (longwave_net  ,scale( 5), 42)
      call read1 (rflvs         ,scale( 6), 52)
      call read1 (toa_swd       ,scale( 7), 62)
      call read1 (cloud_alb     ,scale( 8), 72)
      call read1 (surf_temp     ,scale( 9), 82)
      call read1 (pw            ,scale(10), 92)
      call read1 (rfclr         ,scale(11),102)
      close (2,status='keep')

c      Calculate derived radiation parameters.

      data shortwave_up /62500*-999/,
     &     shortwave_net/62500*-999/,
     &     par_up       /62500*-999/,
     &     par_net      /62500*-999/,
     &     longwave_up  /62500*-999/,
     &     longwave_down/62500*-999/,
     &     net_rad      /62500*-999/

      data eps/0.98/, sigma/5.6697e-8/

c      Short wave.

      do 20 j=1,ny
        do 10 i=1,nx
          if ((shortwave_down(i,j).ne.-999).and.
     &        (surf_alb      (i,j).ne.-999)) then
            shortwave_up (i,j) = shortwave_down(i,j) * surf_alb    (i,j)
            shortwave_net(i,j) = shortwave_down(i,j) - shortwave_up(i,j)
          endif
   10   continue
   20 continue

c      PAR.

      do 40 j=1,ny
        do 30 i=1,nx
          if ((par_down  (i,j).ne.-999).and.
     &        (surf_alb_v(i,j).ne.-999)) then
            par_up (i,j) = par_down(i,j) * surf_alb_v(i,j)
            par_net(i,j) = par_down(i,j) - par_up    (i,j)
          endif
   30   continue
   40 continue

c      Long wave and net radiation.

      do 60 j=1,ny
        do 50 i=1,nx
          if ((surf_alb_v  (i,j).ne.-999).and.
     &        (surf_temp   (i,j).ne.-999).and.
     &        (longwave_net(i,j).ne.-999)) then
            longwave_up  (i,j) = eps*sigma*(surf_temp(i,j) + 273.15)**4
            longwave_down(i,j) = longwave_net (i,j) + longwave_up (i,j)
            net_rad      (i,j) = shortwave_net(i,j) + longwave_net(i,j)
          endif
   50   continue
   60 continue

      write (*,*) 'Date/time ',julianday,hour,min
      i=nx/2
      j=ny/2
      write (*,'(a,f10.3)') 'Latitude        ',lat           (i,j)
      write (*,'(a,f10.3)') 'Longitude       ',lon           (i,j)
      write (*,'(a,f10.3)') 'Short wave down ',shortwave_down(i,j)
      write (*,'(a,f10.3)') 'par_down        ',par_down      (i,j)
      write (*,'(a,f10.3)') 'surf_alb_v      ',surf_alb_v    (i,j)
      write (*,'(a,f10.3)') 'longwave_net    ',longwave_net  (i,j)
      write (*,'(a,f10.3)') 'rflvs           ',rflvs         (i,j)
      write (*,'(a,f10.3)') 'toa_swd         ',toa_swd       (i,j)
      write (*,'(a,f10.3)') 'cloud_alb       ',cloud_alb     (i,j)
      write (*,'(a,f10.3)') 'surf_temp       ',surf_temp     (i,j)
      write (*,'(a,f10.3)') 'pw              ',pw            (i,j)
      write (*,'(a,f10.3)') 'rfclr           ',rfclr         (i,j)
      write (*,'(a,f10.3)') 'shortwave_up    ',shortwave_up  (i,j)
      write (*,'(a,f10.3)') 'shortwave_net   ',shortwave_net (i,j)
      write (*,'(a,f10.3)') 'par_up          ',par_up        (i,j)
      write (*,'(a,f10.3)') 'par_net         ',par_net       (i,j)
      write (*,'(a,f10.3)') 'longwave_up     ',longwave_up   (i,j)
      write (*,'(a,f10.3)') 'longwave_down   ',longwave_down (i,j)
      write (*,'(a,f10.3)') 'net_rad         ',net_rad       (i,j)


          
c--------------USER CODE HERE------------------

      stop

      end


      subroutine header(jday,ihour,minute)

      save
      character page(60)*81

      read (2,rec=1) page
      read (page(23),'(26x,2i2)') ihour,minute
      read (page(24),'(26x,i3 )') jday

c      Correct the minutes value.

      minute=minute+15

      return

      end

        
      subroutine read1(val,scale,irec1)

c      Read and rescale 1 array.  
c      Note that n1*n2 = 6250*10 = 250*250

      implicit none
      save
      integer         n1     ,n2
      parameter      (n1=6250,n2=10)
      real      val  (n1     ,n2   )
      integer*2 ival (n1           )
      real      scale
      integer   i1,i2,irec1,irec2

      irec2=irec1-1
      do 20 i2=1,n2
        irec2=irec2+1
        read (2,rec=irec2) ival
        do 10 i1=1,n1
          if (ival(i1).gt.-29000) then
            val(i1,i2)=ival(i1)*scale
          else
            val(i1,i2)=-999
          endif
   10   continue
   20 continue

      return
      end
