text:
        program ldcmsursubset

	implicit none
        integer(2), allocatable :: band(:,:)
        integer(2), allocatable :: aerob1(:,:)
        integer(2), allocatable :: aerob2(:,:)
        integer(2), allocatable :: aerob3(:,:)
        integer(2), allocatable :: aerob4(:,:)
        integer(2), allocatable :: aerob5(:,:)
        integer(2), allocatable :: aerob6(:,:)
        integer(2), allocatable :: aerob7(:,:)
        integer(2), allocatable :: pband(:,:)
        integer(2), allocatable :: sband(:,:,:)
	integer(2), allocatable :: tband(:,:,:)
	integer(2), allocatable :: opband(:,:)
	integer(2), allocatable :: prband(:,:)
	integer(2), allocatable :: pgband(:,:)
	integer(2), allocatable :: pbband(:,:)
	integer(2), allocatable :: oband(:,:)
	integer(2), allocatable :: tpband(:,:)
	integer(2), allocatable :: aotband(:,:)
	real tmpb10,tmpb11
	real, allocatable :: tlat(:,:)
	real, allocatable :: tlon(:,:)
	real, allocatable :: toz(:,:)
	real, allocatable :: twv(:,:)
	real, allocatable :: twvi(:,:)
	real, allocatable :: tozi(:,:)
	real, allocatable :: txcmg(:,:)
	real, allocatable :: tycmg(:,:)
	real, allocatable :: tp(:,:)
	real, allocatable :: taero(:,:)
	real, allocatable :: taeros(:,:)
	real, allocatable :: taero3(:,:)
	real, allocatable :: tclears(:,:)
	real, allocatable :: taerost(:,:)
	real, allocatable :: teps(:,:)
	real, allocatable :: tepss(:,:)
	real, allocatable :: tresi(:,:)
	real, allocatable :: slprb1(:,:)
	real, allocatable :: intrb1(:,:)
	real, allocatable :: slprb2(:,:)
	real, allocatable :: intrb2(:,:)
	real, allocatable :: slprb7(:,:)
	real, allocatable :: intrb7(:,:)
	integer (2), allocatable :: ttts(:,:)
	integer (2), allocatable :: tttv(:,:)
	integer (2), allocatable :: ttfi(:,:)
	integer (2), allocatable :: tratiob1(:,:)
	integer (2), allocatable :: tratiob2(:,:)
	integer (2), allocatable :: tratiob7(:,:)
	integer (2), allocatable :: tnit(:,:)
	integer(2), allocatable :: wv(:,:)
	BYTE, allocatable :: oz(:,:)
	BYTE, allocatable :: cloud(:,:)
	BYTE, allocatable :: ipflag(:,:)
	BYTE, allocatable :: aerimp(:,:)
	BYTE, allocatable :: smflag(:,:)
	integer(2), allocatable :: dem(:,:)
	integer(2), allocatable :: ratiob1(:,:)
	integer(2), allocatable :: ratiob2(:,:)
	integer(2), allocatable :: ratiob7(:,:)
	integer(2), allocatable :: intratiob1(:,:)
	integer(2), allocatable :: intratiob2(:,:)
	integer(2), allocatable :: intratiob7(:,:)
	integer(2), allocatable :: slpratiob1(:,:)
	integer(2), allocatable :: slpratiob2(:,:)
	integer(2), allocatable :: slpratiob7(:,:)
	integer(2), allocatable :: andwi(:,:)
	integer(2), allocatable :: sndwi(:,:)
	real th1,th2
	real*8 mall
	integer nit
	integer jj

	character(200) filename(12),fname,filenameanc,filenamehdf

	character(2) suffix(12)
	character padding
	character*8 padding1
	character padding2(552)
	integer*2 pad1,pad2
	integer*2 pad3,pad4
	integer offsettiff
	integer sdsind(12)
	data sdsind /0,1,2,3,4,5,6,-1,7,8,9,10/
	data suffix /"01","02","03","04","05","06","07","08","09","10",
     s           "11","QA"/
	integer ii,nr,nc,ib,als,i,j,nrp,ncp,ierr
	integer nrcmg,nccmg,icmg,jcmg
	real u,v
	real xcmg,ycmg
	real xndwi
	real xts,xfs
	real scalefactor,offset
	integer ihdf
	integer iband
! INPUT PARAMETER
         character*6 adate
	 integer imonth,iday,iyear
	 integer sfsattr
! BEGIN OF HDF PARAMETER BLOCK
	 character*4 sdsname
         integer sfstart, sfselect, sfrdata, sfendacc, sfend,set_qamap,set_proj
	 integer sfscatt,sfwdata,sfcreate,sfginfo
         integer sd_id, sd_id2,sds_id, sds_index, status
         integer start(5), edges(5), stride(5)
         integer nstart(2), nedges(2), nstride(2)
         integer DFACC_READ,DFACC_RDWR,DFNT_CHAR8,DFACC_CREATE
	 integer DFNT_INT16,DFNT_FLOAT32,DFNT_UINT8
         parameter (DFACC_READ = 1)
         parameter (DFACC_RDWR = 3)
         parameter (DFACC_CREATE = 4)
         parameter (DFNT_INT16 = 22)
         parameter (DFNT_UINT8 = 21)
          parameter (DFNT_CHAR8 = 4)
        parameter (DFNT_FLOAT32 = 5)
	 character*80 sds_name
	 integer rank,data_type
	 integer n_attrs
	 integer dim_sizes(5)
	 integer dims(2)
         integer dim_length(5), comp_type, comp_prm(4)
! END OF HDF PARAMETER BLOCK

!begin block atmospheric correction variables
c   Aerosol model (string) and length of the string
       character*80 CAMOD
       integer iendarg
C Look up table for atmospheric and geometric quantities
       real tauray(16)
       real oztransa(16)
       real wvtransa(16),wvtransb(16),wvtransc(16)
       real ogtransa0(16),ogtransa1(16)
       real ogtransb0(16),ogtransb1(16)
       real ogtransc0(16),ogtransc1(16)
       real rolutt(16,7,22,8000),nbfi(22,20),tsmax(22,20),tsmin(22,20)
       real ttv(22,20),nbfic(22,20)
       real tts(22)
       integer indts(22)
       real transt(16,7,22,22)
       real sphalbt(16,7,22),normext(16,7,22)
       real xtsmin,xtsstep,xtvmin,xtvstep,pi
       real aot550nm(22)
       data aot550nm /0.01,0.05,0.10,0.15,0.20,0.30,0.40,0.60,
     &                0.80,1.00,1.20,1.40,1.60,1.80,2.00,
     &                2.30,2.60,3.00,3.50,4.00,4.50,5.00/
       real tpres(7)
       data tpres/1050.,1013.0,900.0,800.0,700.0,600.0,500.0/
c usefull variables
       character*6 sbandname(16)
       character*80 err_msg
C The following arguments are all names of the LUTs to look up.
       character*256 tauraynm
       character*256 gscoefnm
       character*256 anglehdf
       character*256 intrefnm
       character*256 transmnm
       character*256 spheranm
        real rotoa,xtv,xfi,raot550nm,uoz,uwv,pres
!	integer ib
        integer retval
	real roslamb
       data (sbandname(i),i=1,8)/"ldcmb1","ldcmb2","ldcmb3","ldcmb4",
     s  "ldcmb5","ldcmb6","ldcmb7","ldcmb8"/
        real tgo,roatm,ttatmg,satm,xrorayp,next
c       integer ldcmind(9)
c       data ldcmind /9,10,4,1,2,6,7,4,6/
!end 	block atmospheric correction variables
!block sharpening
       real r,g,b,xi,xh,xs
!end block sharpening

       real x0,y0,x,y,gsize,lat,lon
       integer row0,col0,utmzone,row,col
       real pres11,pres12,pres21,pres22
       integer uoz11,uoz21,uoz12,uoz22

       real erelc(8),troatm(16)
       real btgo(8),broatm(8),bttatmg(8),bsatm(8)
       real tbtgo(16,22),tbroatm(16,22),tbttatmg(16,22),tbsatm(16,22)
       integer iband1,iband3
       real raot,residual,rsurf
       real fac,xmus,corf,corfb
       integer k,l
       real xcals,xcalo,k1b10,k2b10,k1b11,k2b11,tmpf
       data xcals/3.3420E-04/
       data xcalo/0.10000/
       data k1b10/774.89/
       data k1b11/480.89/
       data k2b10/1321.08/
       data k2b11/1201.14/
       real lat1,lat2,lon1,lon2
       integer colp,rowp
       real dy,dx,ang
       integer nbclear
       real*8 anom,stemp,mclear
       real fack,facl,cldh
       integer cldhmin,cldhmax,icldh
       integer mband5,mband5k,mband5l
       real tcloud
       integer*2 cir,cld,clda,snow,clds,cldt,wat
       real cfac
c bit location of weight for cloudmask QA
       data cld/1/
       data cir/0/
       data clda/2/
       data clds/3/
       data cldt/4/
       data wat/7/
       data cfac/6.0/
       real aaot,sresi
       real fndvi
       integer nbaot,step,hole
       real ros4,ros5
       character*512 qa
       integer rstep,rstepc
       integer iverbose
       integer frln,lrln,frcl,lrcl
       integer fpln,lpln,fpcl,lpcl
       integer nci,nri
	integer xd,yd,zone,sphere
	real uplx,uply,lorx,lory,wbc,ebc,nbc,sbc
	real ros2b1
	integer supind,infind,intst,inten,x1,x2
	real rb1,rb2
	real slpr11,slpr12,slpr21,slpr22
	real intr11,intr12,intr21,intr22
	integer isuccess
	character*200 pfileout,fileout
	character*22 sstr
	integer iout
	integer iaot,ifast,iaots
	real eps
	integer ieps
	real epsres(-1:1)
	real epsaot(-1:1)
	real epsval(-1:1)
	real xa,xb,xc,xdp,xe,xf,coefa,coefb,epsopt
	real taeroavg,tepsavg
	integer nbaeroavg
	real tclearavg
	integer nbclearavg
	real dist,sumdist,smaot
	integer ic,jc,ip1,im1,jp1,jm1
	integer nbpixnf,ipass,iwind
	real xndsi
	real eps1,eps2,eps3,residual1,residual2,residual3,epsmin
	real sraot1,sraot2,sraot3
	real rosband1
	real c1i,c1j,c2i,c2j,c3i,c3j,c4i,c4j,sci,scj,rotang
	integer flag
	real hour,xmin,sec
	integer doy
	real dt,xtu,tu
        integer(2) vza(6366),vaa(6366)
	real saa,sza,viewzenith,viewazimuth
	integer ri,rj
	integer wind,nbpixtot

!initialisation for look up table
	eps=1.
        iaots=1
  	ifast=0
       iverbose=0
       xtv=0.
       xfi=0.
       xtsmin=0
       xtsstep=4.0
       xtvmin=2.84090
       xtvstep=(6.52107-2.84090)
	CAMOD="URBANCLEAN-V2.0"
	iendarg=10
	tauraynm="LDCMLUT/tauray-ldcm.ASC"
	gscoefnm="LDCMLUT/gascoef-ldcm.ASC"
	anglehdf="LDCMLUT/ANGLE_NEW.hdf"
	intrefnm="LDCMLUT/RES_LUT_V3.0-URBANCLEAN-V2.0.hdf"
	transmnm="LDCMLUT/TRANS_LUT_V3.0-URBANCLEAN-V2.0.ASCII"
	spheranm="LDCMLUT/AERO_LUT_V3.0-URBANCLEAN-V2.0.ASCII"
	call readluts(CAMOD,iendarg,tauray,oztransa,wvtransa,
     s                      wvtransb,wvtransc,ogtransa0,ogtransa1,
     s                      ogtransb0,ogtransb1,ogtransc0,ogtransc1,
     s		            tsmax,tsmin,ttv,tts,nbfi,nbfic,indts,
     s                      rolutt,transt,sphalbt,normext,sbandname,err_msg,
     s                      retval,tauraynm,gscoefnm,anglehdf,intrefnm,
     s                      transmnm, spheranm)
        write(6,*) "the luts for urban clean case v2.0 have been read"
	write(6,*) "we can now perform atmospheric correction"

c read the DEM
	nc=7200
	nr=3600
       allocate (dem(nc,nr),stat=ierr)
       sd_id= sfstart("CMGDEM.hdf",DFACC_READ)
       start(1)=0
       start(2) = 0
       edges(1) = nc
       edges(2) = nr
       stride(1) = 1
       stride(2) = 1
       sds_index = 0
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,dem)
       write(6,*) "status", status
       if (status.eq.-1) then
       write(6,*) "Error reading CMGDEM.hdf "
       stop
       endif
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status
c close HDF file
       status = sfend(sd_id)
       write(6,*) "status sfend ",status

       write(6,*) "DEM READ ", dem(2001,1001)

c read the RATIOFILE
	nc=7200
	nr=3600
       allocate (ratiob1(nc,nr),stat=ierr)
       allocate (ratiob2(nc,nr),stat=ierr)
       allocate (ratiob7(nc,nr),stat=ierr)
       allocate (intratiob1(nc,nr),stat=ierr)
       allocate (intratiob2(nc,nr),stat=ierr)
       allocate (intratiob7(nc,nr),stat=ierr)
       allocate (slpratiob1(nc,nr),stat=ierr)
       allocate (slpratiob2(nc,nr),stat=ierr)
       allocate (slpratiob7(nc,nr),stat=ierr)
       allocate (andwi(nc,nr),stat=ierr)
       allocate (sndwi(nc,nr),stat=ierr)
      sd_id= sfstart("ratiomapndwiexp.hdf",DFACC_READ)
      start(1)=0
       start(2) = 0
       edges(1) = nc
       edges(2) = nr
       stride(1) = 1
       stride(2) = 1

       sds_index = 3
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,ratiob2)
       write(6,*) "status", status
       if (status.eq.-1) then
       write(6,*) " Error reading ratiomapndwiexp.hdf"
       stop
       endif

       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

       sds_index = 2
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,ratiob1)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

       sds_index = 4
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,ratiob7)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

       sds_index = 14
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,sndwi)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

       sds_index = 21
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,slpratiob1)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status
       sds_index = 22
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,intratiob1)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

       sds_index = 24
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,slpratiob2)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status
       sds_index = 25
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,intratiob2)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

       sds_index = 27
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,slpratiob7)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status
       sds_index = 28
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,intratiob7)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

       sds_index = 6
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,andwi)
       write(6,*) "status", status
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

c close HDF file
       status = sfend(sd_id)
       write(6,*) "status sfend ",status

       write(6,*) "RATIO READ ", ratiob1(2001,1001),ratiob2(2001,1001),ratiob7(2001,1001)
       write(6,*) "RATIO READ ", slpratiob1(2001,1001),slpratiob2(2001,1001),slpratiob7(2001,1001)
       write(6,*) "RATIO READ ", intratiob1(2001,1001),intratiob2(2001,1001),intratiob7(2001,1001)
       write(6,*) "RATIO READ ", andwi(2001,1001)

cccc GAIL OUT of place
	pi=atan(1.)*4.
	read(5,*) ihdf
	if (ihdf.eq.2) then
	read(5,'(A200)') filenamehdf
	read(5,'(A200)') filenameanc
	read(5,*) xts,xfs
	read(5,*) utmzone,row0,col0,y0,x0
	frln=1
	frcl=1
	write(fileout,'(A200)') "correcteddata.hdf"
c
	else
c the input format is tiff
c read the geometry file
        sd_id= sfstart("l8geom.hdf",DFACC_READ)
        write(6,*) "sd_id",sd_id
	nr=1
	nc=6366
        start(1)=0
        start(2) = 0
        edges(1) = nc
        edges(2) = nr
        stride(1) = 1
        stride(2) = 1
        sds_index = 0
        sds_id    = sfselect(sd_id, sds_index)
        status = sfrdata(sds_id, start, stride, edges,vza)
        write(6,*) "status", status
	if (status.ne.0) then
	write(6,*) "error reading l8geom.hdf"
	stop
	endif
        status=sfendacc(sds_id)
 	write(6,*) "status sfendacc ",status
        sds_index = 1
        sds_id    = sfselect(sd_id, sds_index)
        status = sfrdata(sds_id, start, stride, edges,vaa)
        write(6,*) "status", status
        status=sfendacc(sds_id)
 	write(6,*) "status sfendacc ",status
        status = sfend(sd_id)
c read the input files
	do i=1,12
	read(5,'(A200)') filename(i)
	enddo
	read(5,'(A100)') filenameanc
	read(5,*) nr,nc,nrp,ncp
	read(5,*) xts,xfs
	read(5,*) utmzone,row0,col0,y0,x0
	read(5,*) doy,hour,xmin,sec
	tu=hour+(xmin/60.)+(sec/60./60.)
	read(5,'(A200)') pfileout
	if (ihdf.eq.3) then
	read(5,*) frln,lrln,frcl,lrcl
	read(5,*) c1i,c1j,c2i,c2j,c3i,c3j,c4i,c4j
	read(5,*) offsettiff
	else
	offsettiff=8
	read(5,*,end=19) frln,lrln,frcl,lrcl
	endif
 19     continue

        if (frln.eq.0) then
	write(fileout,'(A200)') pfileout
	frln=1
	lrln=nr
	frcl=1
	lrcl=nc
	iout=0
	else
	ii=index(pfileout," ")-5
	write(sstr,'(A2,4(A1,I4.4))') "ss","-",frln,"-",lrln,"-",frcl,"-",lrcl
	write(fileout,*) pfileout(1:ii),sstr,".hdf"
	iout=2
	endif
	endif

	gsize=30.
	raot550nm=0.06
        fac=pi/180.
	xmus=cos(xts*fac)
cccc END - GAIL OUT of place

C Read ozone and water vapor
	nrcmg=3600
	nccmg=7200
       allocate (oz(nccmg,nrcmg),stat=ierr)
       allocate (wv(nccmg,nrcmg),stat=ierr)
       ii=index(filenameanc," ")-1
        fname=filenameanc(1:ii)
       sd_id= sfstart(fname,DFACC_READ)
      start(1)=0
       start(2) = 0
       edges(1) = nccmg
       edges(2) = nrcmg
       stride(1) = 1
       stride(2) = 1

       sds_index = 0
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,oz)
       write(6,*) "status", status
       if (status.eq.-1) Then
       write(6,*) "Error reading ",fname
       stop
       endif
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

       sds_index = 1
       sds_id    = sfselect(sd_id, sds_index)
       write(6,*) "sds_id", sds_id
       status = sfrdata(sds_id, start, stride, edges,wv)
       write(6,*) "status", status
       if (status.eq.-1) Then
       write(6,*) "Error reading ",fname
       stop
       endif
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status

c close HDF file
       status = sfend(sd_id)
       write(6,*) "status sfend ",status

       write(6,*) "Ozone and Water Vapor read"

c
c	allocate(band(nr,nc),STAT=als)
c	allocate(pband(nrp,ncp),STAT=als)
c	allocate(prband(nrp,ncp),STAT=als)
c	allocate(pgband(nrp,ncp),STAT=als)
c	allocate(pbband(nrp,ncp),STAT=als)
c	allocate(opband(ncp,nrp),STAT=als)
c	allocate(sband(12,nr,nc),STAT=als)
c	allocate(tband(12,nr,nc),STAT=als)
c	allocate(oband(nc,nr),STAT=als)
c	allocate(aotband(nc,nr),STAT=als)
! Getting parameter for atmospheric correction
!        write(6,*) " aot550nm,pressure [Millibars] ,uoz [cm.atm],uwv [g/cm2]"
C update to get the parameter of the scene center
        raot550nm=0.12
	pres=1013.0
	uoz=0.30
	uwv=0.5
	xtv=0.
	xfi=0.
 	write(6,*) raot550nm,pres,uoz,uwv
	write(6,*) xts,xtv,xfi

c open HDF subset
        if  (ihdf.eq.2) then
        ii=index(filenamehdf," ")-1
        fname=filenamehdf(1:ii)
        sd_id= sfstart(fname,DFACC_READ)
        write(6,*) "sd_id",sd_id
        sds_index = 0
        sds_id    = sfselect(sd_id, sds_index)
        write(6,*) "sds_id", sds_id
        status= sfginfo(sds_id, sds_name, rank, dim_sizes, data_type,n_attrs)
       write(6,*) "sdsname ",sds_name
       status = sfendacc(sds_id)
       write(6,*) "status sfendacc ",status
       write(6,*) dim_sizes(1),dim_sizes(2)
       nc= dim_sizes(1)
       nr=dim_sizes(2)
	frln=1
	lrln=nr
	frcl=1
	lrcl=nc
	row0=row0-nr/2
	col0=col0-nc/2
       endif

       if (ihdf.eq.3) then
       nci=lrcl-frcl+1
       nri=lrln-frln+1
       else
       nci=nc
       nri=nr
       endif

       allocate(band(nci,nri),STAT=als)
       write(6,*) "allocate status ",als
c change nc,nr to subset size
       nr=lrln-frln+1
       nc=lrcl-frcl+1
       rstep=nc/50
       allocate(aerob1(nc,nr),STAT=als)
       write(6,*) "allocate status ",als
       allocate(aerob2(nc,nr),STAT=als)
       write(6,*) "allocate status ",als
       allocate(aerob3(nc,nr),STAT=als)
       write(6,*) "allocate status ",als
       allocate(aerob4(nc,nr),STAT=als)
       write(6,*) "allocate status ",als
       allocate(aerob5(nc,nr),STAT=als)
       write(6,*) "allocate status ",als
       allocate(aerob6(nc,nr),STAT=als)
       write(6,*) "allocate status ",als
       allocate(aerob7(nc,nr),STAT=als)
       write(6,*) "allocate status ",als
       allocate(sband(12,nc,nr),STAT=als)
       allocate(tband(12,nc,nr),STAT=als)
       allocate(oband(nc,nr),STAT=als)
       allocate(tpband(nc,nr),STAT=als)
       allocate(aotband(nc,nr),STAT=als)
       allocate(tlat(nc,nr),STAT=als)
       allocate(tlon(nc,nr),STAT=als)
       allocate(txcmg(nc,nr),STAT=als)
       allocate(tycmg(nc,nr),STAT=als)
       allocate(twv(nc,nr),STAT=als)
       allocate(twvi(nc,nr),STAT=als)
       allocate(tozi(nc,nr),STAT=als)
       allocate(tp(nc,nr),STAT=als)
       allocate(taero(nc,nr),STAT=als)
       allocate(taeros(nc,nr),STAT=als)
       allocate(tepss(nc,nr),STAT=als)
       allocate(tclears(nc,nr),STAT=als)
       allocate(taerost(nc,nr),STAT=als)
       allocate(taero3(nc,nr),STAT=als)
       allocate(teps(nc,nr),STAT=als)
       allocate(tresi(nc,nr),STAT=als)
       allocate(tratiob1(nc,nr),STAT=als)
       allocate(tratiob2(nc,nr),STAT=als)
       allocate(tratiob7(nc,nr),STAT=als)
       allocate(tnit(nc,nr),STAT=als)
       allocate(cloud(nc,nr),STAT=als)
       allocate(ipflag(nc,nr),STAT=als)
       allocate(aerimp(nc,nr),STAT=als)
       allocate(smflag(nc,nr),STAT=als)
       allocate(slprb1(nc,nr),STAT=als)
       allocate(intrb1(nc,nr),STAT=als)
       allocate(slprb2(nc,nr),STAT=als)
       allocate(intrb2(nc,nr),STAT=als)
       allocate(slprb7(nc,nr),STAT=als)
       allocate(intrb7(nc,nr),STAT=als)
c saving angles if subset or validation mode (ihdf=3)
      if ((ihdf.eq.3).or.(frln.ne.0)) then
       allocate(ttts(nc,nr),STAT=als)
       allocate(tttv(nc,nr),STAT=als)
       allocate(ttfi(nc,nr),STAT=als)
       endif
c initialize cloud mask
      do i=1,nc
      do j=1,nr
      cloud(i,j)=0
      enddo
      enddo

       if (ihdf.eq.2) then
       start(1)=0
       start(2) = 0
       edges(1) = nci
       edges(2) = nri
       stride(1) = 1
       stride(2) = 1

c read band 12 first
        ib=12
	sds_index=sdsind(ib)
        sds_id =sfselect(sd_id, sds_index)
        write(6,*) "sds_id", sds_id
        status = sfrdata(sds_id, start, stride, edges,band)
         write(6,*) "status", status
         status=sfendacc(sds_id)
         write(6,*) "status sfendacc ",status
	 else
	 ib=12
	 ii=index(filename(ib)," ")-1
	 write(6,*) "reading ",filename(ib)(1:ii)
c compute the offset
          if (ihdf.eq.3) then
	  open(1,file=filename(ib)(1:ii),form='UNFORMATTED',action='READ',access='DIRECT',recl=2*nci*nri+offsettiff)
	  read(1,rec=1) (padding,i=1,offsettiff),((band(i,j),i=1,nci),j=1,nri)
	  close(1)
	  write(6,*) "QA BAND  150,150",band(150,150)
	  else
	  open(1,file=filename(ib)(1:ii),form='UNFORMATTED',action='READ',access='DIRECT',recl=2*nci*nri+offsettiff)
	  read(1,rec=1) (padding,i=1,offsettiff),((band(i,j),i=1,nci),j=1,nri)
	  close(1)
	 endif
	 endif
c
c look for the 4 corners of the image
       if (ihdf.eq.3) goto 561
       if (ihdf.eq.2) goto 561

c read band 5 to determine image rotation
          ib=5
	  open(1,file=filename(ib)(1:ii),form='UNFORMATTED',action='READ',access='DIRECT',recl=2*nci*nri+offsettiff)
	  read(1,rec=1) (padding,i=1,offsettiff),((band(i,j),i=1,nci),j=1,nri)
	  close(1)

       flag=0
       j=1
       i=1
       do while ((flag.eq.0).and.(j.le.nri))
       if (band(i,j).ne.0) then
       write(6,*) "first corner row,col",j,i
       c1j=j
       c1i=i
       flag=1
       endif
       i=i+1
       if (i.gt.nci) then
       i=1
       j=j+1
       endif
       enddo

       flag=0
       j=1
       i=nci
       do while ((flag.eq.0).and.(i.ge.1))
       if (band(i,j).ne.0) then
       write(6,*) "second corner row,col",j,i
       c2j=j
       c2i=i
       flag=1
       endif
       j=j+1
       if (j.gt.nri) then
       j=1
       i=i-1
       endif
       enddo

       flag=0
       j=nri
       i=1
       do while ((flag.eq.0).and.(j.ge.1))
       if (band(i,j).ne.0) then
       write(6,*) "third corner row,col",j,i
       c3j=j
       c3i=i
       flag=1
       endif
       i=i+1
       if (i.gt.nci) then
       i=1
       j=j-1
       endif
       enddo

       flag=0
       j=nri
       i=1
       do while ((flag.eq.0).and.(i.le.nci))
       if (band(i,j).ne.0) then
       write(6,*) "fourth corner row,col",j,i
       c4j=j
       c4i=i
       flag=1
       endif
       j=j-1
       if (j.le.1) then
       j=nri
       i=i+1
       endif
       enddo
c restore band12
          ib=12
	  open(1,file=filename(ib)(1:ii),form='UNFORMATTED',action='READ',access='DIRECT',recl=2*nci*nri+offsettiff)
	  read(1,rec=1) (padding,i=1,offsettiff),((band(i,j),i=1,nci),j=1,nri)
	  close(1)

  561  continue
c
C compute scene center and image rotation
       sci=(c1i+c2i+c3i+c4i)/4.
       scj=(c1j+c2j+c3j+c4j)/4.
       write(6,*) "scene center row,col",scj,sci
       rotang=atan((c2j-c1j)/(c2i-c1i))*180./pi
       write(6,*) "scene rotation degree",rotang

         if (ihdf.eq.3) then
         do i=1,nci
	 do j=1,nri
	 sband(ib,i,j)=band(i,j)
	 enddo
	 enddo
	 else
         do i=frcl,lrcl
	 do j=frln,lrln
	 k=i-frcl+1
	 l=j-frln+1
	 sband(ib,k,l)=band(i,j)
	 enddo
	 enddo
	 endif

c using scene center to compute atmospheric parameter
	 row=frln+nr/2+row0
	 col=frcl+nc/2+col0
c

	 call utmtodeg(utmzone,row,col,x0,y0,gsize,lat,lon)
	 ycmg=(89.975-lat)/0.05+1.
	 xcmg=(179.975+lon)/0.05+1
	 icmg=int(ycmg+0.5)
	 jcmg=int(xcmg+0.5)
         if (wv(jcmg,icmg).ne.0) then
	 uwv=wv(jcmg,icmg)/200.
	 else
	 uwv=0.5
	 endif

         if (oz(jcmg,icmg).ne.0) then
	 uoz=oz(jcmg,icmg)/400.
	 else
	 uoz=0.3
	 endif

         if (dem(jcmg,icmg).ne.-9999) then
	 pres=1013.*exp(-dem(jcmg,icmg)/8500.)
	 else
	 pres=1013.
	 endif
	 raot550nm=0.05

	do ib=1,11
	if (ib.lt.9) then
	iband=ib
         do iaot=1,22
	 raot550nm=aot550nm(iaot)
	 eps=2.5
	 call atmcorlamb2(xts,xtv,xfi,raot550nm,iband,pres,tpres,
     s       aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       uoz,uwv,tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,
     s       rotoa,roslamb,tgo,roatm,ttatmg,satm,xrorayp,next,
     s       err_msg,retval,eps)

	     tbtgo(ib,iaot)=tgo
	     tbroatm(ib,iaot)=roatm
	     tbttatmg(ib,iaot)=ttatmg
	     tbsatm(ib,iaot)=satm
         enddo
	  tgo=tbtgo(ib,2)
	  roatm=tbroatm(ib,2)
	  ttatmg=tbttatmg(ib,2)
	  satm=tbsatm(ib,2)
         endif
!computing parameter for atmospheric correction
	if (ib.ne.8) then
	   if (ihdf.eq.2) then
	   sds_index=sdsind(ib)
           sds_id =sfselect(sd_id, sds_index)
           write(6,*) "sds_id", sds_id
           status = sfrdata(sds_id, start, stride, edges,tpband)
           write(6,*) "status", status
           status=sfendacc(sds_id)
           write(6,*) "status sfendacc ",status
           else
 	  ii=index(filename(ib)," ")-1
	  write(6,*) "reading ",filename(ib)(1:ii)
	  if (ihdf.eq.3) then
	  open(1,file=filename(ib)(1:ii),form='UNFORMATTED',action='READ',access='DIRECT',recl=2*nci*nri+offsettiff)
	  read(1,rec=1) (padding,i=1,offsettiff),((band(i,j),i=1,nci),j=1,nri)
	  close(1)
	  else
	  open(1,file=filename(ib)(1:ii),form='UNFORMATTED',action='READ',access='DIRECT',recl=2*nci*nri+offsettiff)
	  read(1,rec=1) (padding,i=1,offsettiff),((band(i,j),i=1,nci),j=1,nri)
	  close(1)
          endif
c calibrate band 1 to 9 except 8
          if (ihdf.eq.3) then
          if (ib.le.9) then
             do i=1,nci
	     do j=1,nri
 	     if (band(i,j).ge.0) then
             rotoa=(band(i,j)*2.0000E-05)-0.1
	     else
	     rotoa=((65536+band(i,j))*2.0000E-05)-0.1
	     endif
	     row=j+row0-1+frln
	     col=i+col0-1+frcl
	     call utmtodeg(utmzone,row,col,x0,y0,gsize,lat,lon)
             ri=int(6366/2+(col-sci)*cos(-rotang*pi/180.)-(row-scj)*sin(-rotang*pi/180.))
             rj=int((col-sci)*sin(-rotang*pi/180.)+(row-scj)*cos(-rotang*pi/180.))
             dt=-(rj*32./6524.)/60./60.
             xtu=tu+dt
             call possolp(doy,xtu,lon,lat,sza,saa)
             xmus=cos(sza*pi/180.)
	     tpband(i,j)=int(rotoa*10000./xmus)
	     enddo
	     enddo
	     endif
c (ib.eq.9)
	  else
c (ihdf eq 3)
          if (ib.le.9) then
             do i=frcl,lrcl
	     do j=frln,lrln
 	     if (band(i,j).ge.0) then
             rotoa=(band(i,j)*2.0000E-05)-0.1
	     else
	     rotoa=((65536+band(i,j))*2.0000E-05)-0.1
	     endif
     	     k=i-frcl+1
	     l=j-frln+1
             row=j+row0-1
	     col=i+col0-1
	     call utmtodeg(utmzone,row,col,x0,y0,gsize,lat,lon)
             ri=int(6366/2+(col-sci)*cos(-rotang*pi/180.)-(row-scj)*sin(-rotang*pi/180.))
             rj=int((col-sci)*sin(-rotang*pi/180.)+(row-scj)*cos(-rotang*pi/180.))
             dt=-(rj*32./6524.)/60./60.
             xtu=tu+dt
             call possolp(doy,xtu,lon,lat,sza,saa)
             xmus=cos(sza*pi/180.)
	     tpband(k,l)=int(rotoa*10000./xmus)
	     enddo
	     enddo
	     endif
c ib.eq.9
          endif
c ihdf eq 3
	  endif
c ib ne 8

c call the atmospheric correction
             if (ib.lt.9) then
	     	do i=1,nc
	     	do j=1,nr
	     	if (sband(12,i,j).ne.1) then
             	   rotoa=tpband(i,j)/10000.
	           if (ib.eq.1) aerob1(i,j)=tpband(i,j)
	           if (ib.eq.2) aerob2(i,j)=tpband(i,j)
	           if (ib.eq.3) aerob3(i,j)=tpband(i,j)
	           if (ib.eq.4) aerob4(i,j)=tpband(i,j)
	           if (ib.eq.5) aerob5(i,j)=tpband(i,j)
	           if (ib.eq.6) aerob6(i,j)=tpband(i,j)
	           if (ib.eq.7) aerob7(i,j)=tpband(i,j)
                   roslamb=rotoa/tgo
 	           roslamb=roslamb-roatm
 	           roslamb=roslamb/ttatmg
 	           roslamb=roslamb/(1.+satm*roslamb)
	           sband(ib,i,j)=int(roslamb*10000.)
	           btgo(ib)=tgo
	           broatm(ib)=roatm
	           bttatmg(ib)=ttatmg
     	           bsatm(ib)=satm
	       else
	           sband(ib,i,j)=-1000
	       endif
	       enddo
	       enddo
	     else
	        if (ihdf.eq.2) then
                   do i=frcl,lrcl
	           do j=frln,lrln
     	           k=i-frcl+1
	           l=j-frln+1
	           sband(ib,k,l)=tpband(i,j)
	           enddo
	           enddo
	        else
		   if (ihdf.eq.3) then
                      do i=1,nc
	              do j=1,nr
  	              sband(ib,i,j)=band(i,j)
	              enddo
	              enddo
		   else
                      do i=frcl,lrcl
	              do j=frln,lrln
     	              k=i-frcl+1
	              l=j-frln+1
	              sband(ib,k,l)=band(i,j)
	              enddo
	              enddo
		   endif
	         endif
	   endif

	   if (ib.eq.9) then
	   do i=1,nc
	   do j=1,nr
	   sband(ib,i,j)=tpband(i,j)
	   enddo
	   enddo
	   endif

	else
c iband eq 8
	write(6,*) "not reading band8 sorry"
	endif

c
	enddo

c calibrate thermal bands
         do i=1,nc
	 do j=1,nr
	 if (sband(10,i,j).lt.0)  then
	    tmpb10=65536+sband(10,i,j)
	    else
	    tmpb10=sband(10,i,j)
	    endif

	 if (sband(11,i,j).lt.0)  then
	     tmpb11=65536+sband(11,i,j)
	     else
	     tmpb11=sband(11,i,j)
	     endif

	 if ((sband(12,i,j).ne.1)) then
	 tmpf=xcals*tmpb10+xcalo
	 tmpf=k2b10/log(1.+k1b10/tmpf)
	 sband(10,i,j)=int(tmpf*10.)
	 tmpf=xcals*tmpb11+xcalo
	 tmpf=k2b11/log(1.+k1b11/tmpf)
	 sband(11,i,j)=int(tmpf*10.)
	 endif
	 enddo
	 enddo
         write(6,*) "Finished calibrating Thermal"

c compute lat,lon,rowcmg,colcmg,ozone,water vapor,pressure
c skipping reading the aerosol product goto 111
         do i=1,nc
c	 write(6,*) " Processing column ",i
	 if ((i-int(i/rstep)*rstep).eq.0) then
	 write(6,*) " processing collumn ",i
	 write(6,*) " iaots ",iaots
	 endif
	 do j=1,nr
c	 goto 11
	 if ((sband(12,i,j).eq.1)) then
	 ipflag(i,j)=5
	 goto 11
	 endif
	 if ((j.eq.673).and.(i.eq.709)) then
	 write(6,*) "J'en tiens un"
	 endif

	 row=frln+(j-1)+row0
	 col=frcl+(i-1)+col0
	 call utmtodeg(utmzone,row,col,x0,y0,gsize,lat,lon)
	 tlat(i,j)=lat
	 tlon(i,j)=lon
	 ycmg=(89.975-lat)/0.05+1.
	 xcmg=(179.975+lon)/0.05+1
	 txcmg(i,j)=xcmg
	 tycmg(i,j)=ycmg
c	 icmg=int(ycmg+0.5)
c	 jcmg=int(xcmg+0.5)
c         twv(i,j)=wv(jcmg,icmg)
	 icmg=int(ycmg)
	 jcmg=int(xcmg)
	 u=(ycmg-icmg)
	 v=(xcmg-jcmg)
	 twvi(i,j)=wv(jcmg,icmg)*(1.-u)*(1.-v)+wv(jcmg+1,icmg)*(1.-u)*v
     s	      +wv(jcmg,icmg+1)*(1.-v)*u+wv(jcmg+1,icmg+1)*u*v
         twvi(i,j)=twvi(i,j)/100.

	 if (oz(jcmg,icmg).lt.0) then
	 uoz11=256+oz(jcmg,icmg)
	 else
	 uoz11=oz(jcmg,icmg)
	 endif

	 if (oz(jcmg+1,icmg).lt.0) then
	 uoz12=256+oz(jcmg+1,icmg)
	 else
	 uoz12=oz(jcmg+1,icmg)
	 endif

	 if (oz(jcmg,icmg+1).lt.0) then
	 uoz21=256+oz(jcmg,icmg+1)
	 else
	 uoz21=oz(jcmg,icmg+1)
	 endif

	 if (oz(jcmg+1,icmg+1).lt.0) then
	 uoz22=256+oz(jcmg+1,icmg+1)
	 else
	 uoz22=oz(jcmg+1,icmg+1)
	 endif

	 if (uoz11.eq.0) uoz11=120
	 if (uoz12.eq.0) uoz12=120
	 if (uoz21.eq.0) uoz21=120
	 if (uoz22.eq.0) uoz22=120

	 tozi(i,j)=uoz11*(1.-u)*(1.-v)+uoz12*(1.-u)*v
     s	      +uoz21*(1.-v)*u+uoz22*u*v

         tozi(i,j)=tozi(i,j)/400.

	 if (dem(jcmg,icmg).ne.-9999) then
	 pres11=1013.*exp(-dem(jcmg,icmg)/8500.)
	 else
	 pres11=1013.
	 cloud(i,j)=-128
c	 tresi(i,j)=-1.
	 endif
	 if (dem(jcmg+1,icmg).ne.-9999) then
	 pres12=1013.*exp(-dem(jcmg+1,icmg)/8500.)
	 else
	 pres12=1013.
	 endif
	 if (dem(jcmg,icmg+1).ne.-9999) then
	 pres21=1013.*exp(-dem(jcmg,icmg+1)/8500.)
	 else
	 pres21=1013.
	 endif
	 if (dem(jcmg+1,icmg+1).ne.-9999) then
	 pres22=1013.*exp(-dem(jcmg+1,icmg+1)/8500.)
	 else
	 pres22=1013.
	 endif
	 tp(i,j)=pres11*(1.-u)*(1.-v)+pres12*(1.-u)*v
     s	      +pres21*(1.-v)*u+pres22*u*v

c inverting aerosol
	do ib=1,8
	erelc(ib)=-1.
	enddo
	rb1=ratiob1(jcmg,icmg)/1000.
	rb2=ratiob2(jcmg,icmg)/1000.
	if ((rb2.gt.1.0).or.(rb1.gt.1.0).or.(rb2.lt.0.1).or.(rb1.lt.0.1)) then
	slpratiob1(jcmg,icmg)=0
	slpratiob2(jcmg,icmg)=0
	slpratiob7(jcmg,icmg)=0
	intratiob1(jcmg,icmg)=550
	intratiob2(jcmg,icmg)=600
	intratiob7(jcmg,icmg)=2000
	else
	if (sndwi(jcmg,icmg).lt.200) then
	slpratiob1(jcmg,icmg)=0
	slpratiob2(jcmg,icmg)=0
	slpratiob7(jcmg,icmg)=0
	intratiob1(jcmg,icmg)=ratiob1(jcmg,icmg)
	intratiob2(jcmg,icmg)=ratiob2(jcmg,icmg)
	intratiob7(jcmg,icmg)=ratiob7(jcmg,icmg)
	endif
	endif

	rb1=ratiob1(jcmg+1,icmg)/1000.
	rb2=ratiob2(jcmg+1,icmg)/1000.
	if ((rb2.gt.1.0).or.(rb1.gt.1.0).or.(rb2.lt.0.1).or.(rb1.lt.0.1)) then
	slpratiob1(jcmg+1,icmg)=0
	slpratiob2(jcmg+1,icmg)=0
	slpratiob7(jcmg+1,icmg)=0
	intratiob1(jcmg+1,icmg)=550
	intratiob2(jcmg+1,icmg)=600
	intratiob7(jcmg+1,icmg)=2000
	else
	if (sndwi(jcmg+1,icmg).lt.200) then
	slpratiob1(jcmg+1,icmg)=0
	slpratiob2(jcmg+1,icmg)=0
	slpratiob7(jcmg+1,icmg)=0
	intratiob1(jcmg+1,icmg)=ratiob1(jcmg+1,icmg)
	intratiob2(jcmg+1,icmg)=ratiob2(jcmg+1,icmg)
	intratiob7(jcmg+1,icmg)=ratiob7(jcmg+1,icmg)
	endif
	endif

	rb1=ratiob1(jcmg,icmg+1)/1000.
	rb2=ratiob2(jcmg,icmg+1)/1000.
	if ((rb2.gt.1.0).or.(rb1.gt.1.0).or.(rb2.lt.0.1).or.(rb1.lt.0.1)) then
	slpratiob1(jcmg,icmg+1)=0
	slpratiob2(jcmg,icmg+1)=0
	slpratiob7(jcmg,icmg+1)=0
	intratiob1(jcmg,icmg+1)=550
	intratiob2(jcmg,icmg+1)=600
	intratiob7(jcmg,icmg+1)=2000
	else
	if (sndwi(jcmg,icmg+1).lt.200) then
	slpratiob1(jcmg,icmg+1)=0
	slpratiob2(jcmg,icmg+1)=0
	slpratiob7(jcmg,icmg+1)=0
	intratiob1(jcmg,icmg+1)=ratiob1(jcmg,icmg+1)
	intratiob2(jcmg,icmg+1)=ratiob2(jcmg,icmg+1)
	intratiob7(jcmg,icmg+1)=ratiob7(jcmg,icmg+1)
	endif
	endif

	rb1=ratiob1(jcmg+1,icmg+1)/1000.
	rb2=ratiob2(jcmg+1,icmg+1)/1000.

	if ((rb2.gt.1.0).or.(rb1.gt.1.0).or.(rb2.lt.0.1).or.(rb1.lt.0.1)) then
	slpratiob1(jcmg+1,icmg+1)=0
	slpratiob2(jcmg+1,icmg+1)=0
	slpratiob7(jcmg+1,icmg+1)=0
	intratiob1(jcmg+1,icmg+1)=550
	intratiob2(jcmg+1,icmg+1)=600
	intratiob7(jcmg+1,icmg+1)=2000
	else
	if (sndwi(jcmg+1,icmg+1).lt.200) then
	slpratiob1(jcmg+1,icmg+1)=0
	slpratiob2(jcmg+1,icmg+1)=0
	slpratiob7(jcmg+1,icmg+1)=0
	intratiob1(jcmg+1,icmg+1)=ratiob1(jcmg+1,icmg+1)
	intratiob2(jcmg+1,icmg+1)=ratiob2(jcmg+1,icmg+1)
	intratiob7(jcmg+1,icmg+1)=ratiob7(jcmg+1,icmg+1)
	endif
	endif

	 slpr11=slpratiob2(jcmg,icmg)/1000.
	 intr11=intratiob2(jcmg,icmg)/1000.
	 slpr12=slpratiob2(jcmg+1,icmg)/1000.
	 intr12=intratiob2(jcmg+1,icmg)/1000.
	 slpr21=slpratiob2(jcmg,icmg+1)/1000.
	 intr21=intratiob2(jcmg,icmg+1)/1000.
	 slpr22=slpratiob2(jcmg+1,icmg+1)/1000.
	 intr22=intratiob2(jcmg+1,icmg+1)/1000.
	 slprb2(i,j)=slpr11*(1.-u)*(1.-v)+slpr12*(1.-u)*v
     s	      +slpr21*(1.-v)*u+slpr22*u*v
	 intrb2(i,j)=intr11*(1.-u)*(1.-v)+intr12*(1.-u)*v
     s	      +intr21*(1.-v)*u+intr22*u*v

	 slpr11=slpratiob1(jcmg,icmg)/1000.
	 intr11=intratiob1(jcmg,icmg)/1000.
	 slpr12=slpratiob1(jcmg+1,icmg)/1000.
	 intr12=intratiob1(jcmg+1,icmg)/1000.
	 slpr21=slpratiob1(jcmg,icmg+1)/1000.
	 intr21=intratiob1(jcmg,icmg+1)/1000.
	 slpr22=slpratiob1(jcmg+1,icmg+1)/1000.
	 intr22=intratiob1(jcmg+1,icmg+1)/1000.
	 slprb1(i,j)=slpr11*(1.-u)*(1.-v)+slpr12*(1.-u)*v
     s	      +slpr21*(1.-v)*u+slpr22*u*v
	 intrb1(i,j)=intr11*(1.-u)*(1.-v)+intr12*(1.-u)*v
     s	      +intr21*(1.-v)*u+intr22*u*v

	 slpr11=slpratiob7(jcmg,icmg)/1000.
	 intr11=intratiob7(jcmg,icmg)/1000.
	 slpr12=slpratiob7(jcmg+1,icmg)/1000.
	 intr12=intratiob7(jcmg+1,icmg)/1000.
	 slpr21=slpratiob7(jcmg,icmg+1)/1000.
	 intr21=intratiob7(jcmg,icmg+1)/1000.
	 slpr22=slpratiob7(jcmg+1,icmg+1)/1000.
	 intr22=intratiob7(jcmg+1,icmg+1)/1000.
	 slprb7(i,j)=slpr11*(1.-u)*(1.-v)+slpr12*(1.-u)*v
     s	      +slpr21*(1.-v)*u+slpr22*u*v
	 intrb7(i,j)=intr11*(1.-u)*(1.-v)+intr12*(1.-u)*v
     s	      +intr21*(1.-v)*u+intr22*u*v

        xndwi=dble(sband(5,i,j)-dble(sband(7,i,j)/2.))/dble(sband(5,i,j)+dble(sband(7,i,j)/2.))
	th1=((andwi(jcmg,icmg)+2.*sndwi(jcmg,icmg))/1000.)
	th2=((andwi(jcmg,icmg)-2.*sndwi(jcmg,icmg))/1000.)
	if (xndwi.gt.th1) xndwi=th1
	if (xndwi.lt.th2) xndwi=th2
        erelc(1)=(xndwi*dble(slprb1(i,j))+intrb1(i,j))
        erelc(2)=(xndwi*dble(slprb2(i,j))+intrb2(i,j))
        erelc(7)=(xndwi*dble(slprb7(i,j))+intrb7(i,j))
 	erelc(4)=1.
	tratiob1(i,j)=int(erelc(1)*1000.)
	tratiob2(i,j)=int(erelc(2)*1000.)
	tratiob7(i,j)=int(erelc(7)*1000.)

	troatm(1)=aerob1(i,j)/10000.
 	troatm(2)=aerob2(i,j)/10000.
	troatm(4)=aerob4(i,j)/10000.
	troatm(7)=aerob7(i,j)/10000.
	iband1=4
	iband3=1
	pres=tp(i,j)
	uoz=tozi(i,j)
	uwv=twvi(i,j)
	if (((i.eq.2).and.(j.eq.1)).or.((i.eq.10000).and.(j.eq.10000))) then
	write(6,*) "aeroband1 ",aerob1(i,j)
       write(6,*) i,j,xts,xtv,xfi,pres,uoz,uwv,erelc(1),erelc(2),erelc(7),
     s	troatm(1),troatm(2),troatm(4),troatm(7),iband1,iband3,tratiob1(i,j)
        iverbose=1
	else
	iverbose=0
        endif

       ifast=0
       eps=1.0
       iaots=1
	 row=frln+(j-1)+row0
	 col=frcl+(i-1)+col0
	 call utmtodeg(utmzone,row,col,x0,y0,gsize,lat,lon)
       ri=int(6366/2+(col-sci)*cos(-rotang*pi/180.)-(row-scj)*sin(-rotang*pi/180.))
       rj=int((col-sci)*sin(-rotang*pi/180.)+(row-scj)*cos(-rotang*pi/180.))
       dt=-(rj*32./6524.)/60./60.
       if (ri.lt.1) ri=1
       if (ri.gt.6366) ri=6366
       viewzenith=vza(ri)
       viewazimuth=vaa(ri)
       xtu=tu+dt
       call possolp(doy,xtu,lon,lat,sza,saa)
       xtv=viewzenith/100.
       xts=sza
       xfi=saa-(viewazimuth/100.)
       if ((ihdf.eq.3).or.(frln.ne.0)) then
       ttts(i,j)=int(xts*100.)
       tttv(i,j)=int(xtv*100.)
       ttfi(i,j)=int(xfi*100.)
       endif
       call subaeroretv3(iband1,iband3,xts,xtv,xfi,pres,uoz,uwv,erelc,troatm,
     c       tpres,aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,iverbose,raot,residual,nit,next,ierr,
     s       ifast,tbtgo,tbroatm,tbttatmg,tbsatm,iaots,eps)
     	 if (iverbose.eq.1) then
         write(6,*) "aero retrieved i,j,raot,eps,residual" ,i,j,raot,eps,residual
	 endif
	 eps1=eps
	 residual1=residual
	 sraot1=raot

       eps=1.75
       call subaeroretv3(iband1,iband3,xts,xtv,xfi,pres,uoz,uwv,erelc,troatm,
     c       tpres,aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,iverbose,raot,residual,nit,next,ierr,
     s       ifast,tbtgo,tbroatm,tbttatmg,tbsatm,iaots,eps)
     	 if (iverbose.eq.1) then
         write(6,*) "aero retrieved i,j,raot,eps,residual" ,i,j,raot,eps,residual
	 endif
	 eps2=eps
	 residual2=residual
	 sraot2=raot

       eps=2.5
       call subaeroretv3(iband1,iband3,xts,xtv,xfi,pres,uoz,uwv,erelc,troatm,
     c       tpres,aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,iverbose,raot,residual,nit,next,ierr,
     s       ifast,tbtgo,tbroatm,tbttatmg,tbsatm,iaots,eps)
     	 if (iverbose.eq.1) then
         write(6,*) "aero retrieved i,j,raot,eps,residual" ,i,j,raot,eps,residual
	 endif
	 eps3=eps
	 residual3=residual
	 sraot3=raot

c find eps that minimize the residual
	xa=(eps1*eps1)-(eps3*eps3)
	xd=(eps2*eps2)-(eps3*eps3)
    	xb=(eps1-eps3)
	xe=(eps2-eps3)
	xc=residual1-residual3
	xf=residual2-residual3
	coefa=(xc*xe-xb*xf)/(xa*xe-xb*xd)
	coefb=(xa*xf-xc*xd)/(xa*xe-xb*xd)
	epsmin=-coefb/(2*coefa)
       eps=epsmin
       if ((epsmin.ge.1.0).and.(epsmin.le.2.5)) then
       eps=epsmin
       call subaeroretv3(iband1,iband3,xts,xtv,xfi,pres,uoz,uwv,erelc,troatm,
     c       tpres,aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,iverbose,raot,residual,nit,next,ierr,
     s       ifast,tbtgo,tbroatm,tbttatmg,tbsatm,iaots,eps)
         else
	 if (epsmin.le.1.0) then
	 eps=eps1
	 residual=residual1
	 raot=sraot1
	 endif
	 if (epsmin.ge.2.5) then
	 eps=eps3
	 residual=residual3
	 raot=sraot3
	 endif
	 endif

     	 if (iverbose.eq.1) then
         write(6,*) "aero retrieved i,j,raot,eps,residual" ,i,j,raot,eps,residual
	 endif

c	 endif
         if (ierr.eq.1) then
	 write(6,*) "problem in subaeroretv3 at i,j ",i,j
	 endif

         tnit(i,j)=nit
	 taero3(i,j)=next*raot
	 teps(i,j)=eps
         corf=raot/xmus
	 taero(i,j)=raot

c         if (residual.lt.((0.015+0.005*corf)*aerob7(i,j)/300.)) then
c         if (troatm(7).gt.0.05) then
c	 corfb=troatm(7)/0.05
c	 else
	 corfb=1.
c	 endif
         if (residual.lt.((0.015+0.005*corf+0.10*troatm(7))*corfb)) then
c test if band5 makes sense

	iband=5
	rotoa=aerob5(i,j)/10000.
	raot550nm=raot
	call atmcorlamb2(xts,xtv,xfi,raot550nm,iband,pres,tpres,
     s       aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       uoz,uwv,tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,
     s       rotoa,roslamb,tgo,roatm,ttatmg,satm,xrorayp,next,
     s       err_msg,retval,eps)
        ros5=roslamb
	iband=4
	rotoa=aerob4(i,j)/10000.
	raot550nm=raot
	call atmcorlamb2(xts,xtv,xfi,raot550nm,iband,pres,tpres,
     s       aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       uoz,uwv,tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,
     s       rotoa,roslamb,tgo,roatm,ttatmg,satm,xrorayp,next,
     s       err_msg,retval,eps)
         ros4=roslamb
         if ((ros5.gt.0.1).and.((ros5-ros4)/(ros4+ros5).gt.0)) then
	 taero(i,j)=raot
	 tresi(i,j)=residual
	 ipflag(i,j)=0
	 else
c this is water the retrieval needs to be redone
	 taero(i,j)=raot

	 tresi(i,j)=residual
	 ipflag(i,j)=-1
	 endif

	 else
	 taero(i,j)=raot
	 tresi(i,j)=residual
	 ipflag(i,j)=-1
	 endif

c redo the retrieval if water
        if (ipflag(i,j).eq.-1) then
	do ib=1,8
	erelc(ib)=-1.
	enddo
	troatm(1)=aerob1(i,j)/10000.
	troatm(4)=aerob4(i,j)/10000.
	troatm(5)=aerob5(i,j)/10000.
	troatm(7)=aerob7(i,j)/10000.
	iband1=4
        erelc(7)=1.0
        erelc(5)=1.0
        erelc(4)=1.0
        erelc(1)=1.0

	pres=tp(i,j)
	uoz=tozi(i,j)
	uwv=twvi(i,j)
	if (((i.eq.408).and.(j.eq.254))) then
	write(6,*) "aeroband1 ",aerob1(i,j)
       write(6,*) i,j,xts,xtv,xfi,pres,uoz,uwv,erelc(1),erelc(2),erelc(7),
     s	troatm(1),troatm(2),troatm(4),troatm(7),iband1,iband3,tratiob1(i,j)
        iverbose=1
	else
	iverbose=0
        endif

       ifast=0

       eps=1.5
       iaots=1
       call subaeroretwat(iband1,iband3,xts,xtv,xfi,pres,uoz,uwv,erelc,troatm,
     c       tpres,aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,iverbose,raot,residual,nit,next,ierr,
     s       ifast,tbtgo,tbroatm,tbttatmg,tbsatm,iaots,eps)
     	 if (iverbose.eq.1) then
         write(6,*) "aero retrieved i,j,raot,eps,residual" ,i,j,raot,eps,residual
	 endif

     	 if (iverbose.eq.1) then
         write(6,*) "aero retrieved i,j,raot,eps,residual" ,i,j,raot,eps,residual
	 endif

c	 endif
         tnit(i,j)=nit
	 taero3(i,j)=next*raot
	 teps(i,j)=eps
         corf=raot/xmus
	 taero(i,j)=raot
	 tresi(i,j)=residual
c test band 1 reflectance eliminate negative
	iband=1
	rotoa=aerob1(i,j)/10000.
	raot550nm=raot
	call atmcorlamb2(xts,xtv,xfi,raot550nm,iband,pres,tpres,
     s       aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       uoz,uwv,tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,
     s       rotoa,roslamb,tgo,roatm,ttatmg,satm,xrorayp,next,
     s       err_msg,retval,eps)
        rosband1=roslamb

         if ((residual.gt.((0.010+0.005*corf))).or.(rosband1.lt.0)) then
c 	 write(6,*) "problem in subaeroretwat at i,j ",i,j
         ipflag(i,j)=2
          else
         tnit(i,j)=nit
	 taero3(i,j)=next*raot
	 teps(i,j)=eps
         corf=raot/xmus
	 taero(i,j)=raot
	 tresi(i,j)=residual
	 endif
  	 endif

 11      continue
	 enddo
	 enddo

CCCCCC GAIL HERE CCCC
c before interpolation increase the interpolation to neighbourring pixels (3x3?)
c not sure this is needed
c       goto 1001
       do j=1,nr
       do i=1,nc
       if ((ipflag(i,j).gt.1).and.(ipflag(i,j).lt.5)) then
       do k=i-2,i+2
       do l=j-2,j+2
       if ((k.ge.1).and.(l.ge.1).and.(k.le.nc).and.(l.le.nr)) then
       if (ipflag(k,l).eq.0) ipflag(k,l)=9
       endif
       enddo
       enddo
       endif
       enddo
       enddo

       do j=1,nr
       do i=1,nc
       if (ipflag(i,j).eq.9) ipflag(i,j)=2
       enddo
       enddo
c
c average the aerosol 10x10 for valid value
       nbpixnf=0
       nbpixtot=0
       do j=1,nr
       do i=1,nc
       smflag(i,j)=0
       if(sband(12,i,j).ne.1) then
       nbpixtot=nbpixtot+1
       taeroavg=0.
       tepsavg=0.
       nbaeroavg=0
       do k=max(j-5,1),min(j+5,nr)
       do l=max(i-5,1),min(i+5,nc)
c       if (j.le.5) write(6,*) "i,j,k ",i,j,k
       if ((ipflag(l,k).eq.0).or.(ipflag(l,k).eq.-1)) then
       nbaeroavg=nbaeroavg+1
       taeroavg=taero(l,k)+taeroavg
       tepsavg=teps(l,k)+tepsavg
       endif
       enddo
       enddo
c       if (j.le.5) write(6,*) "i,j,nbaeroavg",i,j,nbaeroavg
       if (nbaeroavg.gt.20) then
       taeroavg=taeroavg/nbaeroavg
       tepsavg=tepsavg/nbaeroavg
       taeros(i,j)=taeroavg
       tepss(i,j)=tepsavg
       smflag(i,j)=1
       else
c       write(6,*) "Pixels at i,j not filled",i,j
       nbpixnf=nbpixnf+1
       endif
       endif
       enddo
       enddo

       if (nbpixnf.eq.nbpixtot) then
        do j=1,nr
       do i=1,nc
       if(sband(12,i,j).ne.1) then
       taeros(i,j)=0.05
       tepss(i,j)=1.5
       smflag(i,j)=1
       endif
       enddo
       enddo
       goto 121
       endif

       write(6,*) "Second Pass"
c while nbpixnf.ne.0
       ipass=2
       wind=5
       do while (nbpixnf.ne.0)
       write(6,*) " Pass number ",ipass, "nbpixnf ",nbpixnf
       nbpixnf=0
       do j=1,nr
       do i=1,nc
       if((sband(12,i,j).ne.1).and.(smflag(i,j).eq.0)) then
       tepsavg=0.
       taeroavg=0.
       nbaeroavg=0
       do k=max(j-wind,1),min(j+wind,nr)
       do l=max(i-wind,1),min(i+wind,nc)
       if (smflag(l,k).eq.1) then
       nbaeroavg=nbaeroavg+1
       taeroavg=taeros(l,k)+taeroavg
       tepsavg=tepss(l,k)+tepsavg
      endif
       enddo
       enddo
       if (nbaeroavg.gt.0) then
       taeroavg=taeroavg/nbaeroavg
       tepsavg=tepsavg/nbaeroavg
       taeros(i,j)=taeroavg
       tepss(i,j)=tepsavg
       smflag(i,j)=1
       else
c       write(6,*) "Second pass Pixels at i,j not filled",i,j
       nbpixnf=nbpixnf+1
       endif
       endif
       enddo
       enddo
       ipass=ipass+1
c       wind=wind+10
       enddo

c end of averaging

c
 121   continue
       do j=1,nr
       do i=1,nc
       if ((i.eq.408).and.(j.eq.254)) then
       write(6,*) "verbose on in interpolation ",i,j,ipflag(i,j)
       iverbose=1
       else
       iverbose=0
       endif

        if ((ipflag(i,j).gt.1).and.(.not.btest(cloud(i,j),cir)).and.(.not.btest(cloud(i,j),cld))
     s    .and.(ipflag(i,j).lt.5)) then
       taero(i,j)=taeros(i,j)
       teps(i,j)=tepss(i,j)
       ipflag(i,j)=1
       endif
       enddo
       enddo

 2005  continue
c

c perform the correction
 1001     continue
	  rstepc=rstep*5
          do ib=1,7
	  do i=1,nc
	 if ((i-int(i/rstepc)*rstepc).eq.0) then
	 write(6,*) " atmospheric correction collumn ",i," band ",ib
	 endif
	  do j=1,nr
	 row=frln+(j-1)+row0
	 col=frcl+(i-1)+col0
	 call utmtodeg(utmzone,row,col,x0,y0,gsize,lat,lon)
       ri=int(6366/2+(col-sci)*cos(-rotang*pi/180.)-(row-scj)*sin(-rotang*pi/180.))
       rj=int((col-sci)*sin(-rotang*pi/180.)+(row-scj)*cos(-rotang*pi/180.))
       dt=-(rj*32./6524.)/60./60.
       if (ri.lt.1) ri=1
       if (ri.gt.6366) ri=6366
       viewzenith=vza(ri)
       viewazimuth=vaa(ri)
       xtu=tu+dt
       call possolp(doy,xtu,lon,lat,sza,saa)
       xtv=viewzenith/100.
       xts=sza
       xfi=saa-(viewazimuth/100.)

c	  write(6,*) "i , j",i,j,tresi(i,j)
	  if (sband(12,i,j).ne.1) then
	  if ((.not.btest(cloud(i,j),cir)).and.(.not.btest(cloud(i,j),cld))) then
	  rsurf=sband(ib,i,j)/10000.
	  rotoa=(rsurf*bttatmg(ib)/(1.-bsatm(ib)*rsurf)+broatm(ib))*btgo(ib)
	  if (ib.eq.1) then
c	  write(6,*) "band 1", aerob1(i,j),rotoa
	  endif
	  raot550nm=taero(i,j)
	  eps=teps(i,j)
	  pres=tp(i,j)
	  uwv=twvi(i,j)
	  uoz=tozi(i,j)
	  iband=ib
c	  write(6,*) xts,xtv,xfi,raot550nm,iband,pres,uoz,uwv,tauray
	call atmcorlamb2(xts,xtv,xfi,raot550nm,iband,pres,tpres,
     s       aot550nm,rolutt,
     s       transt,xtsstep,xtsmin,xtvstep,xtvmin,
     s       sphalbt,normext,
     s       tsmax,tsmin,nbfic,nbfi,tts,indts,ttv,
     s       uoz,uwv,tauray,
     s       ogtransa0,ogtransa1,ogtransb0,ogtransb1,
     s       ogtransc0,ogtransc1,
     s       wvtransa,wvtransb,wvtransc,oztransa,
     s       rotoa,roslamb,tgo,roatm,ttatmg,satm,xrorayp,next,
     s       err_msg,retval,eps)
          if ((i.eq.676).and.(j.eq.383)) then
	  iverbose=1
	  write(6,*) "  test on negative i,j,rotoa,roslamb ",i,j,rotoa,roslamb," tauaero ",taero(i,j)
	  else
	  iverbose=0
	  endif
          if ((i.eq.2111).and.(j.eq.4123)) then
	  write(6,*) "roslamb ",roslamb," tauaero ",taero(i,j)
	  endif
          if (ib.eq.1) then
c endif or refining optical depth
c set up aerosol QA
c          write(6,*) "SETTING UP aerosol qa BITS and aerosol impact
          aerimp(i,j)=min(int(abs(rsurf-roslamb)*1000),255)
          if (abs(rsurf-roslamb).le.0.015) then
	    cloud(i,j)=cloud(i,j)+16
	    else
	    if   (abs(rsurf-roslamb).lt.0.03) then
	          cloud(i,j)=cloud(i,j)+32
	          else
	          cloud(i,j)=cloud(i,j)+48
	          endif
	    endif
	  endif
C endif of ib=1

	  sband(ib,i,j)=int(roslamb*10000.)
	  endif
	  endif
	  enddo
	  enddo
	  enddo
c update the cloud mask
	   nbclear=0
	   mclear=0.
	   mall=0.

 	  do i=1,nc
 	  do j=1,nr
	  xndsi=(sband(3,i,j)-sband(6,i,j))*1.0/(sband(3,i,j)+sband(6,i,j))
          if (sband(12,i,j).ne.1) then
	  if ((tresi(i,j).gt.0.05).and.((sband(10,i,j)/10.).lt.305.)) then
	      cloud(i,j)=cloud(i,j)+2
	      else
	   nbclear=nbclear+1
	   mclear=mclear+(sband(10,i,j)/10.)
	   endif
	   endif
	  enddo
	  enddo
	   if (nbclear.gt.0) then
	   mclear=mclear/nbclear
	   else
	   mclear=305.0
	   endif
           write(6,*) "average clear temperature  %clear", mclear,nbclear*100./(nr*nc)

c compute temperature of clear pixels
       if (nbclear.eq.0) then
       write(6,*) "nbclear eq 0 "
       do j=1,nr
       do i=1,nc
       tclears(i,j)=mclear
       enddo
       enddo
       goto 123
       endif
c set up adjacent cloud bit
       do i=1,nc
       do j=1,nr
       if (sband(12,i,j).ne.1) then
       if (btest(cloud(i,j),cld).or.(btest(cloud(i,j),cir))) then
       do k=max(i-5,1),min(i+5,nc)
       do l=max(j-5,1),min(j+5,nr)
       if ((k.ge.1).and.(k.le.nc).and.(l.ge.1).and.(l.le.nr)) then
       if ((.NOT.btest(cloud(k,l),cld)).and.(.not.btest(cloud(k,l),cir)).and.(.not.btest(cloud(k,l),clda))) then
          cloud(k,l)=cloud(k,l)+4
	  endif
	endif
	enddo
	enddo
	endif
	endif
	enddo
	enddo

       nbpixnf=0
       do j=1,nr
       do i=1,nc
       smflag(i,j)=0
       if(sband(12,i,j).ne.1) then
       tclearavg=0.
       nbclearavg=0
       do k=max(j-10,1),min(j+10,nr)
       do l=max(i-10,1),min(i+10,nc)
       if ((.not.btest(cloud(l,k),cld)).and.(.not.btest(cloud(l,k),clda)).and.(sband(12,l,k).ne.1)) then
       nbclearavg=nbclearavg+1
       tclearavg=(sband(10,l,k)/10.)+tclearavg
       endif
       enddo
       enddo
       if (nbclearavg.gt.0) then
       tclearavg=tclearavg/nbclearavg
       tclears(i,j)=tclearavg
       smflag(i,j)=1
       else
       nbpixnf=nbpixnf+1
       endif
       endif
       enddo
       enddo

       write(6,*) "Second Pass temperature"
       do j=1,nr
       do i=1,nc
       if(sband(12,i,j).ne.1) then
       if (smflag(i,j).eq.0) then
       tclears(i,j)=mclear
       endif
       endif
       enddo
       enddo
       goto 123

       ipass=2
       iwind=20
       do while (nbpixnf.ne.0)
       write(6,*) " Pass number ",ipass, "nbpixnf ",nbpixnf
       nbpixnf=0
       do j=1,nr
       do i=1,nc
       if((sband(12,i,j).ne.1).and.(smflag(i,j).eq.0)) then
       tclearavg=0.
       nbclearavg=0
       do k=max(j-iwind,1),min(j+iwind,nr)
       do l=max(i-iwind,1),min(i+iwind,nc)
       if (smflag(l,k).eq.1) then
       nbclearavg=nbclearavg+1
       tclearavg=tclears(l,k)+tclearavg
       endif
       enddo
       enddo
       if (nbclearavg.gt.0) then
       tclearavg=tclearavg/nbclearavg
       tclears(i,j)=tclearavg
       smflag(i,j)=1
       else
       nbpixnf=nbpixnf+1
       endif
       endif
       enddo
       enddo
       ipass=ipass+1
       enddo
c refine cloud mask
 123   continue
       do j=1,nr
       do i=1,nc
       if (btest(cloud(i,j),cld)) then
       if (((sband(10,i,j)/10.).ge.(tclears(i,j)-0.5)).or.((sband(6,i,j)+sband(7,i,j)).lt.1000)) then
       cloud(i,j)=cloud(i,j)-2
c test for snow
       xndsi=(sband(3,i,j)-sband(6,i,j))*1.0/(sband(3,i,j)+sband(6,i,j))
       if (xndsi.gt.0.) then
       cloud(i,j)=cloud(i,j)+64
       endif
       endif
       endif
       enddo
       enddo
c do the cloud adjacency
       do i=1,nc
       do j=1,nr
       if (btest(cloud(i,j),clda)) then
          cloud(i,j)=cloud(i,j)-4
	  endif
       enddo
       enddo
       do i=1,nc
       do j=1,nr
       if (btest(cloud(i,j),cld).or.(btest(cloud(i,j),cir))) then
       do k=i-5,i+5
       do l=j-5,j+5
       if ((k.ge.1).and.(k.le.nc).and.(l.ge.1).and.(l.le.nr)) then
       if ((.NOT.btest(cloud(k,l),cld)).and.(.not.btest(cloud(k,l),cir)).and.(.not.btest(cloud(k,l),clda))) then
          cloud(k,l)=cloud(k,l)+4
	  endif
	endif
	enddo
	enddo
	endif
	enddo
	enddo

c cloud shadow
c compute adjustement to the North
c scene center lat,lon
	 row=frln+(nr-1)/2.
	 col=frcl+(nc-1)/2.
	 write(6,*) "row col ",row,col
	 call utmtodeg(utmzone,row,col,x0,y0,gsize,lat1,lon1)
c move 100 pixels north
         write(6,*) "LAT1 LON1",lat1,lon1
         rowp=row-100
	 call utmtodeg(utmzone,rowp,col,x0,y0,gsize,lat2,lon2)
         write(6,*) "LAT2 LON2",lat2,lon2
	 lon=lon1
	 lat=lat2
	 call degtoutm(utmzone,lat,lon,x0,y0,gsize,colp,rowp)
	 write(6,*) "colp, rowp",colp,rowp
	 dy=row-rowp
	 dx=colp-col
	 ang=atan(dx/dy)*180./pi
	 write(6,*) "Adjustement to true North ", ang
c  compute cloud shadow
       write(6,*) "looking for cloud shadow"
       fack=cos(xfs*fac)*tan(xts*fac)/gsize
       facl=sin(xfs*fac)*tan(xts*fac)/gsize
       do i=1,nc
       do j=1,nr
       if (btest(cloud(i,j),cld).or.btest(cloud(i,j),cir)) then
           tcloud=sband(10,i,j)/10.
           cldh=(tclears(i,j)-tcloud)*1000./cfac
	   if (cldh.lt.0.) cldh=0.
	   cldhmin=cldh-1000.
	   cldhmax=cldh+1000.
	   mband5=9999
	   if (cldhmin.lt.0) cldhmin=0.
           do icldh=cldhmin/10,cldhmax/10
	    cldh=icldh*10.
            k=j+fack*cldh
            l=i-facl*cldh
            if ((k.ge.1).and.(k.le.nr).and.(l.ge.1)
     &	    .and.(l.le.nc)) then

            if ((sband(6,l,k).lt.800).and.((sband(3,l,k)-sband(4,l,k)).lt.100)) then
            if (btest(cloud(l,k),cld).or.btest(cloud(l,k),cir)
     &       .or.btest(cloud(l,k),clds)) then
             continue
             else
c store the value of band6 as well as the l and k value
             if (sband(6,l,k).lt.mband5) then
	     mband5=sband(6,l,k)
	     mband5k=k
	     mband5l=l
	     endif
             endif
             endif
             endif
	     enddo
	     if (mband5.lt.9999) then
	     l=mband5l
	     k=mband5k
             cloud(l,k)=cloud(l,k)+(2**clds)
	     endif
         endif
	 enddo
	 enddo

 966     continue

!Saving in an HDF file
         write(6,*) "saving corrected data"
	 ii=index(fileout," ")-1
	 write(6,*) "fileout ",fileout
         sd_id= sfstart(adjustl(fileout),DFACC_CREATE)
	 write(6,*) " filout ",sd_id
	 do ib=1,12
	 if (ib.ne.8) then
	 do i=1,nc
	 do j=1,nr
	 oband(i,j)=sband(ib,i,j)
	 enddo
	 enddo
	 if (ib.eq.2) then
         write(sds_name,'(A4,A2,A5)') "band",suffix(ib),"-blue"
	 endif
	 if (ib.eq.3) then
         write(sds_name,'(A4,A2,A6)') "band",suffix(ib),"-green"
	 endif
	 if (ib.eq.4) then
         write(sds_name,'(A4,A2,A4)') "band",suffix(ib),"-red"
	 endif
	 if ((ib.le.1).or.(ib.ge.5)) then
         write(sds_name,'(A4,A2)') "band",suffix(ib)
	 endif
	 write(6,*) "writing band ",sds_name
	 dim_length(1)=nc
	 dim_length(2)=nr
	 comp_type=4
	 comp_prm(1)=8
	 rank=2
	 dim_sizes(1)=nc
	 dim_sizes(2)=nr
         start(1) = 0
         start(2) = 0
          edges(1) = nc
         edges(2) = nr
         stride(1) = 1
         stride(2) = 1
	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,oband)
c check if QA band to write the attribute
	 write(6,*) "status sfwdata ",status
          if ((ib.ge.1).and.(ib.le.9)) then
	 write(6,*) "writing attribute band ",sds_id
	 scalefactor=10000.
	 offset=0.0
	 status=sfsattr(sds_id, "scale_factor", DFNT_FLOAT32, 1, scalefactor)
	 status=sfsattr(sds_id, "offset", DFNT_FLOAT32, 1, offset)
	endif
          if ((ib.eq.10).or.(ib.eq.11)) then
	 write(6,*) "writing attribute band ",sds_id
	 scalefactor=10.
	 offset=0.0
	 status=sfsattr(sds_id, "scale_factor", DFNT_FLOAT32, 1, scalefactor)
	 status=sfsattr(sds_id, "offset", DFNT_FLOAT32, 1, offset)
	endif
        status = sfendacc(sds_id)
	 write(6,*) "status sfendacc ",status
	 endif
	 enddo
c writing lat and lon
c            write(sds_name,'(A3)') "lat"
c   	 sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
c            status=sfwdata(sds_id,start,stride,edges,tlat)
c   	 write(6,*) "status sfwdata ",status
c            status = sfendacc(sds_id)
c            write(sds_name,'(A3)') "lon"
c   	 sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
c            status=sfwdata(sds_id,start,stride,edges,tlon)
c   	 write(6,*) "status sfwdata ",status
c            status = sfendacc(sds_id)
c writing xcmg and ycmg
c          write(sds_name,'(A4)') "xcmg"
c 	 sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
c          status=sfwdata(sds_id,start,stride,edges,txcmg)
c 	 write(6,*) "status sfwdata ",status
c          status = sfendacc(sds_id)
c          write(sds_name,'(A4)') "ycmg"
c 	 sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
c          status=sfwdata(sds_id,start,stride,edges,tycmg)
c 	 write(6,*) "status sfwdata ",status
c          status = sfendacc(sds_id)
c writing twv,twvi

c        write(sds_name,'(A3)') "twv"
c	 sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
cc       status=sfwdata(sds_id,start,stride,edges,twv)
c	 write(6,*) "status sfwdata ",status
cc       status = sfendacc(sds_id)
         if (ihdf.ne.0) iout=1
         if (iout.ne.0) then
c writing aerosol algorithm bands (toa)
         if (ihdf.eq.3) then
         write(sds_name,'(A7,A2)') "toaband",suffix(1)
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,aerob1)
         status = sfendacc(sds_id)
         write(sds_name,'(A7,A2)') "toaband",suffix(2)
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,aerob2)
         status = sfendacc(sds_id)
         write(sds_name,'(A7,A2)') "toaband",suffix(3)
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,aerob3)
         status = sfendacc(sds_id)
         write(sds_name,'(A7,A2)') "toaband",suffix(4)
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,aerob4)
         status = sfendacc(sds_id)
         write(sds_name,'(A7,A2)') "toaband",suffix(5)
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,aerob5)
         status = sfendacc(sds_id)
         write(sds_name,'(A7,A2)') "toaband",suffix(6)
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,aerob6)
         status = sfendacc(sds_id)
         write(sds_name,'(A7,A2)') "toaband",suffix(7)
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,aerob7)
         status = sfendacc(sds_id)
	 endif

         write(sds_name,'(A4)') "twvi"
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,twvi)
 	 write(6,*) "status sfwdata ",status
         status = sfendacc(sds_id)
         write(sds_name,'(A4)') "tozi"
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,tozi)
 	 write(6,*) "status sfwdata ",status
         status = sfendacc(sds_id)
         write(sds_name,'(A6)') "tpresi"
 	 sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,tp)
 	 write(6,*) "status sfwdata ",status
         status = sfendacc(sds_id)
	 endif
         if  (iout.ne.0) then
	 write(sds_name,'(A11)') "taero-band1"
	sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,taero3)
	write(6,*) "status sfwdata ",status
	 status = sfendacc(sds_id)
	 write(sds_name,'(A5)') "taero"
	sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,taero)
	write(6,*) "status sfwdata ",status
	 status = sfendacc(sds_id)
	 write(sds_name,'(A6)') "taeros"
	sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,taeros)
	write(6,*) "status sfwdata ",status
	 status = sfendacc(sds_id)
	 write(sds_name,'(A7)') "tclears"
	sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,tclears)
	write(6,*) "status sfwdata ",status
	 status = sfendacc(sds_id)
	 write(sds_name,'(A5)') "tresi"
	sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,tresi)
	write(6,*) "status sfwdata ",status
	 status = sfendacc(sds_id)
	 write(sds_name,'(A5)') "tepsa"
	sds_id=sfcreate(sd_id,sds_name,DFNT_FLOAT32,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,teps)
	write(6,*) "status sfwdata ",status
	 status = sfendacc(sds_id)
  	   write(sds_name,'(A8)') "tratiob1"
	   write(6,*) "tratiob1 ",tratiob1(672,709),tratiob1(709,672)
	sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,tratiob1)
	 write(6,*) "status sfwdata ",status
c	 status = sfendacc(sds_id)
	 scalefactor=1000.
	 offset=0.0
	 status=sfsattr(sds_id, "scale_factor", DFNT_FLOAT32, 1, scalefactor)
	 status=sfsattr(sds_id, "offset", DFNT_FLOAT32, 1, offset)
	 status = sfendacc(sds_id)
	   write(sds_name,'(A8)') "tratiob2"
	sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,tratiob2)
	 write(6,*) "status sfwdata ",status
c	 status = sfendacc(sds_id)
	 scalefactor=1000.
	 offset=0.0
	 status=sfsattr(sds_id, "scale_factor", DFNT_FLOAT32, 1, scalefactor)
	 status=sfsattr(sds_id, "offset", DFNT_FLOAT32, 1, offset)
	 status = sfendacc(sds_id)
	   write(sds_name,'(A8)') "tratiob7"
	sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,tratiob7)
	 write(6,*) "status sfwdata ",status
c	 status = sfendacc(sds_id)
	 scalefactor=1000.
	 offset=0.0
	 status=sfsattr(sds_id, "scale_factor", DFNT_FLOAT32, 1, scalefactor)
	 status=sfsattr(sds_id, "offset", DFNT_FLOAT32, 1, offset)
	 status = sfendacc(sds_id)

	 write(sds_name,'(A6)') "solzen"
	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,ttts)
	 write(6,*) "status sfwdata ",status
	 scalefactor=100.
	 offset=0.0
	 status=sfsattr(sds_id, "scale_factor", DFNT_FLOAT32, 1, scalefactor)
	 status=sfsattr(sds_id, "offset", DFNT_FLOAT32, 1, offset)
	 status = sfendacc(sds_id)

	 write(sds_name,'(A7)') "viewzen"
	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,tttv)
	 write(6,*) "status sfwdata ",status
	 scalefactor=100.
	 offset=0.0
	 status=sfsattr(sds_id, "scale_factor", DFNT_FLOAT32, 1, scalefactor)
	 status=sfsattr(sds_id, "offset", DFNT_FLOAT32, 1, offset)
	 status = sfendacc(sds_id)

	 write(sds_name,'(A7)') "relazim"
	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,ttfi)
	 write(6,*) "status sfwdata ",status
	 scalefactor=100.
	 offset=0.0
	 status=sfsattr(sds_id, "scale_factor", DFNT_FLOAT32, 1, scalefactor)
	 status=sfsattr(sds_id, "offset", DFNT_FLOAT32, 1, offset)
	 status = sfendacc(sds_id)

	 endif

	 if (1.eq.2) then
	 write(sds_name,'(A8)') "tratiob1"
	 write(6,*) "tratiob1 ",tratiob1(672,709),tratiob1(709,672)
	 sds_id=sfcreate(sd_id,sds_name,DFNT_INT16,rank,dim_sizes)
	 status=sfwdata(sds_id,start,stride,edges,tratiob1)
	 write(6,*) "status sfwdata ",status
c	 status = sfendacc(sds_id)
	 scalefactor=1000.
	 offset=0.0
	 status=sfsattr(sds_id, "scale_factor", DFNT_FLOAT32, 1, scalefactor)
	 status=sfsattr(sds_id, "offset", DFNT_FLOAT32, 1, offset)
	 status = sfendacc(sds_id)
	 endif

         write(sds_name,'(A5)') "CLOUD"
	 sds_id=sfcreate(sd_id,sds_name,DFNT_UINT8,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,cloud)
	 write(6,*) "status sfwdata ",status
	 write(6,*) "writing attribute data for CLOUD band ",sds_id
	 status=set_qamap(sds_id)
	 write(6,*) "status set qa attribute ",status
         status = sfendacc(sds_id)
         write(sds_name,'(A6)') "IPFLAG"
	 sds_id=sfcreate(sd_id,sds_name,DFNT_UINT8,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,ipflag)
	 write(6,*) "status sfwdata ",status
          write(sds_name,'(A6)') "AERIMP"
	 sds_id=sfcreate(sd_id,sds_name,DFNT_UINT8,rank,dim_sizes)
         status=sfwdata(sds_id,start,stride,edges,aerimp)
	 write(6,*) "status sfwdata ",status
        status = sfendacc(sds_id)
c writing projection metadata
	 xd=nc
	 yd=nr
	 uplx=x0+gsize*(frcl+col0-2)
	 uply=y0-gsize*(frln+row0-2)
	 lorx=uplx+nc*gsize
	 lory=uply-nr*gsize
	 zone=utmzone
	 sphere=12

	 status=set_proj(sd_id,xd,yd,uplx,uply,lorx,lory,zone,sphere,wbc,ebc,nbc,sbc)
	 write(6,*) "status set projection attribute ",status

!closing the HDF file
         status = sfend(sd_id)
	 write(6,*) "status sfend ",status

	stop
	end


text:
	SUBROUTINE GDTXRN ( tbegin, tend, itftyp, ntimef, filtim,
     +			      ntime, times, iret )
C************************************************************************
C* GDTXRN 								*
C*									*
C* This subroutine checks which times in a file are within a time	*
C* range without an increment.						*
C*									*
C* GDTXRN ( TBEGIN, TEND, ITFTYP, NTIMEF, FILTIM, NTIME, TIMES,	IRET )	*
C*									*
C* Input parameters:							*
C*	TBEGIN		CHAR*		Start time			*
C*	TEND		CHAR*		Stop time			*
C*	ITFTYP		INTEGER		Range type			*
C*					  1 = forecast range		*
C*					  2 = date/time range		*
C*					  3 = combo			*
C*	NTIMEF		INTEGER		Number of times in file		*
C*	FILTIM (NTIMEF)	CHAR*		File times			*
C*									*
C* Output parameters:							*
C*	NTIME		INTEGER		Number of output times		*
C*	TIMES (NTIME)	CHAR*		Times in range			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 5/89						*
C************************************************************************
	CHARACTER*(*)	tbegin, tend, filtim (*), times (*)
C------------------------------------------------------------------------
	iret  = 0
	ntime = 0
C
C*	Check times in file depending on ITFTYP.
C
	IF ( itftyp .eq. 1 ) THEN
C
C*	    Get forecast start and end times.
C
	    CALL TG_IFTM ( tbegin (12:12), tbegin (13: ), ifbeg, ier )
	    CALL TG_IFTM ( tend   (12:12), tend   (13: ), ifend, ier )
C
C*	    Check against all times in the file.
C
	    DO i = 1, ntimef
		IF ( filtim (i)( :12 ) .eq. tbegin ( :12 ) )  THEN
		    CALL TG_IFTM ( filtim (i)(12:12), filtim (i)(13: ),
     +				   ifcur, ier )
		    IF ( ( ifcur .ge. ifbeg ) .and.
     +			 ( ifcur .le. ifend ) ) THEN
			ntime = ntime + 1
			times ( ntime ) = filtim (i)
		    END IF
		END IF
	    END DO
C
C*	    Check the case where the date/time is being incremented.
C
	ELSE
	    DO i = 1, ntimef
		IF ( filtim (i)( 12: ) .eq. tbegin ( 12: ) )  THEN
		    IF ( ( filtim (i)( :11 ) .ge. tbegin ( :11 ) ) .and.
     +			 ( filtim (i)( :11 ) .le. tend ( :11 ) ) ) THEN
			ntime = ntime + 1
			times ( ntime ) = filtim (i)
		    END IF
		END IF
	    END DO
	END IF
C*
	IF ( itftyp .eq. 3 ) THEN
	    CALL TG_IFTM ( tend (12:12), tbegin (13:), ifbeg, ier)
	    CALL TG_IFTM ( tend (12:12), tend(13:),    ifend, ier)
	    DO i = 1, ntimef
		IF ( filtim (i)(:12) .eq. tend (:12) ) THEN
		    CALL TG_IFTM ( filtim (i)(12:12), filtim (i)(13:),
     +				   ifcur, ier )
		    IF ( (ifcur .gt. ifbeg) .and.
     +			 (ifcur .le. ifend) ) THEN
			ntime = ntime + 1
			times (ntime) = filtim (i)
		    END IF
		END IF
	    END DO
	END IF
C
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
	SUBROUTINE HSSPLN  ( msltyp, mslstr, msldir, tslsiz, mslwid,
     +			     iret )
C************************************************************************
C* HSSPLN - VG								*
C*									*
C* This subroutine sets the special line attributes.			*
C*									*
C* HSSPLN  ( MSLTYP, MSLSTR, MSLDIR, TSLSIZ, MSLWID, IRET )		*
C*									*
C* Input parameters:							*
C*	MSLTYP		INTEGER		Special line type		*
C*	MSLSTR		INTEGER		Special line stroke multiplier	*
C*	MSLDIR		INTEGER		Special line direction indicator*
C*	TSLSIZ		REAL		Special line size		*
C*	MSLWID		INTEGER		Special line width		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* D. Keiser/GSC	 4/97						*
C* S. Maxwell/GSC        6/97   Documentation changes                   *
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DEVCHR.CMN'
	INCLUDE		'DVWNDW.CMN'
C------------------------------------------------------------------------
	iret = NORMAL
C
	CALL VSSPLN  ( msltyp, mslstr, msldir, tslsiz, mslwid, iret )
C*
	RETURN
	END


text:
module hecmw_pair_array
  use hecmw_util

  private

  public :: hecmwST_pair_array
  public :: hecmw_pair_array_init
  public :: hecmw_pair_array_finalize
  public :: hecmw_pair_array_append
  public :: hecmw_pair_array_sort
  public :: hecmw_pair_array_find_id

  type hecmwST_pair
    integer(kind=kint) :: id
    integer(kind=kint) :: i1, i2
  end type hecmwST_pair

  type hecmwST_pair_array
    integer(kind=kint) :: num
    integer(kind=kint) :: max_num
    type (hecmwST_pair), pointer :: pairs(:) => null()
  end type hecmwST_pair_array

contains

  subroutine hecmw_pair_array_init(parray, max_num)
    implicit none
    type (hecmwST_pair_array), intent(inout) :: parray
    integer(kind=kint), intent(in) :: max_num
    !if (associated(parray%pairs)) deallocate(parray%pairs)
    allocate(parray%pairs(max_num))
    parray%max_num = max_num
    parray%num = 0
  end subroutine hecmw_pair_array_init

  subroutine hecmw_pair_array_finalize(parray)
    implicit none
    type (hecmwST_pair_array), intent(inout) :: parray
    if (associated(parray%pairs)) deallocate(parray%pairs)
    parray%max_num = 0
    parray%num = 0
  end subroutine hecmw_pair_array_finalize

  subroutine hecmw_pair_array_append(parray, id, i1, i2)
    implicit none
    type (hecmwST_pair_array), intent(inout) :: parray
    integer(kind=kint), intent(in) :: id, i1, i2
    if (parray%num >= parray%max_num) then
      stop 'ERROR: hecmw_pair_array_append: overflow'
    endif
    parray%num = parray%num + 1
    parray%pairs(parray%num)%id = id
    parray%pairs(parray%num)%i1 = i1
    parray%pairs(parray%num)%i2 = i2
  end subroutine hecmw_pair_array_append

  subroutine hecmw_pair_array_sort(parray)
    implicit none
    type (hecmwST_pair_array), intent(inout) :: parray
    call pairs_sort(parray%pairs, 1, parray%num)
  end subroutine hecmw_pair_array_sort

  function hecmw_pair_array_find_id(parray, i1, i2)
    implicit none
    integer(kind=kint) :: hecmw_pair_array_find_id
    type (hecmwST_pair_array), intent(inout) :: parray
    integer(kind=kint), intent(in) :: i1, i2
    type (hecmwST_pair) :: p
    integer(kind=kint) :: id
    p%i1 = i1
    p%i2 = i2
    call pairs_find(parray%pairs, 1, parray%num, p, id)
    hecmw_pair_array_find_id = id
  end function hecmw_pair_array_find_id

  function pairs_comp(p1, p2)
    implicit none
    integer(kind=kint) :: pairs_comp
    type (hecmwST_pair), intent(in) :: p1, p2
    if (p1%i1 < p2%i1) then
      pairs_comp = -1
    else if (p1%i1 > p2%i1) then
      pairs_comp = 1
    else
      if (p1%i2 < p2%i2) then
        pairs_comp = -1
      else if (p1%i2 > p2%i2) then
        pairs_comp = 1
      else
        pairs_comp = 0
      endif
    endif
  end function pairs_comp

  recursive subroutine pairs_sort(pairs, from, to)
    implicit none
    type (hecmwST_pair), pointer :: pairs(:)
    integer(kind=kint), intent(in) :: from, to
    integer(kind=kint) :: center, left, right
    type (hecmwST_pair) :: pivot, tmp
    if (from >= to) return
    center = (from + to) / 2
    pivot = pairs(center)
    left = from
    right = to
    do
      do while (pairs_comp(pairs(left), pivot) < 0)
        left = left + 1
      enddo
      do while (pairs_comp(pivot, pairs(right)) < 0)
        right = right - 1
      enddo
      if (left >= right) exit
      tmp = pairs(left)
      pairs(left) = pairs(right)
      pairs(right) = tmp
      left = left + 1
      right = right - 1
    enddo
    if (from < left-1) call pairs_sort(pairs, from, left-1)
    if (right+1 < to) call pairs_sort(pairs, right+1, to)
    return
  end subroutine pairs_sort

  recursive subroutine pairs_find(pairs, from, to, p, id)
    implicit none
    type (hecmwST_pair), pointer :: pairs(:)
    integer(kind=kint), intent(in) :: from, to
    type (hecmwST_pair), intent(in) :: p
    integer(kind=kint), intent(out) :: id
    integer(kind=kint) :: center, icomp
    if (from > to) then
      id = -1
      return
    endif
    center = (from + to) / 2
    icomp = pairs_comp(p, pairs(center))
    if (icomp < 0) then
      call pairs_find(pairs, from, center-1, p, id)
      return
    else if (icomp > 0) then
      call pairs_find(pairs, center+1, to, p, id)
      return
    else
      id = pairs(center)%id
      return
    endif
  end subroutine pairs_find

end module hecmw_pair_array

max_stars_repo_name: fujitsu/FrontISTR

max_issues_repo_name: fujitsu/FrontISTR

max_forks_repo_name: fujitsu/FrontISTR