NASA Ocean Color

ocssw V2022

       subroutine readin(ilm,irh,isd,itau)
 c
 c     subroutine readin reads the opens files and reads input from      
 c     unit 21                                                           
 c***********************************************************************
 c.....include common blocks                                             
       implicit real*8 (a-h,o-z)
       include 'common_all.cmn'
 c***********************************************************************
       integer*4 izx(2)
       real*8 xzx(3),athcnd(25),aphcnd(46),ftrx(16,25,46),
      1       ptix(8,25,46)
 c                                                                       
 c     character*1 ciwind,ciconc
       character*1 ciconc
       character*2 ciwind
       character*2 citau,crh,cset
       character*2 cisd,cilmd,cwind
       character*4 cilm
       character*90 xname1,xname2,xname3,xname4,xname5,file_name
       character*90 xname4a,xname5a,xname6
       character*80 txt,temp
 c                                                                       
       len1=index(dir1,' ')-1
       len2=index(dir2,' ')-1
       len3=index(dir3,' ')-1
       len4=index(dir4,' ')-1
       len5=index(dir5,' ')-1
 c
 c      if(wwl(ilm).lt.1.0)then
          iwl=wwl(ilm)*1.0d3+0.01d0
 c      else
 c         iwl=wwl(ilm)*1.0d4+0.1d0
 c      endif
       call convtc(idnint(wwl(ilm)*1.0d3+0.01),4,cilm)
       write(*,*) wwl(ilm),cilm
       call convtc(isd,2,cisd)
       call convtc(itau,2,citau)
       call convtc(iwind-1,2,ciwind)
 c     call convtc(iconc,1,ciconc)
       call convtc(iww(iwind),2,cwind)
       call convtc(krhum(irh),2,crh)
       call convtc(iset,2,cset)    !iset is passed through the common block zeroa
 c
       if(iref.eq.1)ciwind='00'
       if(iref.eq.0)ciwind='99'
       if(iset.eq.0)crh='00'  !iset is read in sub. inusr and passed through common zeroa
 c
       xname1='rt2_wl'//cilm//'x'//ciwind
       xname2='sd'//cisd//'rh'//crh//'ta'//citau//'_set'//cset
 c     xname3=cilm//ciwind//ciconc//cisd//citau
       xname4a='rt1_'//'wl'//cilm//'sd'//cisd//'rh'//crh//'ta'//citau
       nm4a=index(xname4a,' ')-1
       xname4=xname4a(1:nm4a)//'_set'//cset
       xname5a='rt1_sp'//'wl'//cilm//'sd'//cisd//'rh'//crh//'ta'//citau
       nm5a=index(xname5a,' ')-1
       xname5=xname5a(1:nm5a)//'_set'//cset
 c
 c
       nm1=index(xname1,' ')-1
       nm2=index(xname2,' ')-1
 c     nm3=index(xname3,' ')-1
       nm4=index(xname4,' ')-1
       nm5=index(xname5,' ')-1
 c
 c
       num_arg = iargc()
       if (num_arg .lt. 1) then
         file_name = 'dumpza'
       else
         call getarg(1, file_name)
       end if
 c
       open(3,file=
      1     dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'dn.dat',
      2     status='unknown')
       open(4,file=
      1     dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'up.dat',
      2     status='unknown')
       open(6,file=file_name,status='unknown')
       open(21,file=dir1(1:len1)//'/'//xname4(1:nm4)//'.dat',
      1      status='old')
       open(24,file=
      1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'.dat',
      2 access='direct',recl=1915*4,
      3 form='unformatted',status='unknown')
 c
       open(53,access='direct',
      1     form='unformatted',status='scratch',recl=18400*4) 
       open(54,access='direct',
      1     form='unformatted',status='scratch',recl=18400*4)
       open(55,access='direct',
      1     form='unformatted',status='scratch',recl=18400*4)
       open(64,access='direct',
      1     form='unformatted',status='scratch',recl=18400*4) 
       open(65,access='direct',
      1     form='unformatted',status='scratch',recl=18400*4) 
       open(71,access='direct',
      1     form='unformatted',status='scratch',recl=18400*4) 
       open(72,access='direct',
      1     form='unformatted',status='scratch',recl=18400*4) 
 c
       if(icrft.eq.1)then
       open(13,file=
      1     dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'crftdn.dat',
      2     status='unknown')
       open(14,file=
      1     dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'crftup.dat',
      2     status='unknown')
        open(15,file=
      1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'crft.dat',
      2 access='direct',recl=1915*4,
      3 form='unformatted',status='unknown')
       open(73,access='direct',
      1     form='unformatted',status='scratch',recl=18400*4)
       open(74,access='direct',
      1     form='unformatted',status='scratch',recl=18400*4)
       endif
 c
       if(iactflx.eq.1)then
       open(30,file=
      1     dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'flx.dat',
      2     status='unknown')
       endif
 c
       if(isurf.eq.1)then
       open(16,file=
      1     dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'surfup.dat',
      2     status='unknown')
         if(iref.eq.2)then
            open(41,file=
      1     dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'dirup.dat',
      2     status='unknown')    
            open(42,file=
      1     dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'ocnup.dat',
      2     status='unknown')
            open(43,file=
      1     dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'skyup.dat',
      2     status='unknown')
         endif
       endif
 c                                                                       
 c      rewind(7)
 c      read(7,*)xday,xmonth,xyear
        xday=-777.0
        xmonth=-777.0
        xyear=-777.0
 c
        read(21,125)temp
        read(21,126)wvlth,psrfc,rho,xozn,tautot,deltau,nolyr
        read(21,125)temp
        read(21,126)tr,tm,ta,tcar,twat,tozn,ifc
        read(21,125)temp
        ht(1)=84.0d0
        pl(1)=5.3040d-6
        do i=1,nolyr
           read(21,128)ht(i+1),pl(i+1),dtrr(i),dtmm(i),dtaa(i),dtot(i)
 c         write(6,128)ht(i+1),pl(i+1),dtrr(i),dtmm(i),dtaa(i),dtot(i)
        enddo
        press=pl(nolyr+1)
 c
        if(ipsudo.eq.1)then
             open(2,file=
      1           dir1(1:len1)//'/'//xname5(1:nm5)//'.dat',
      2           status='old')        
             read(2,125)temp
             read(2,145)nmodl
             read(2,125)temp
             do i=1,nmodl                         
                read(2,150)htp(i),ppo(i),taur(i),taum(i),tauabs(i)
             enddo 
 c
       endif
 c                                                                       
       qs=1.0d-15
       qt=1.0d-15
       do 70 i=1,1801
       do 60 j=1,4
          t(i,j) = 0.0d0
    60 continue
    70 continue
       if(ifc.eq.1)then
 c
          xname6='phs_twl'//cilm//'sd'//cisd//'rh'//crh//'_set'//cset
          nm6=index(xname6,' ')-1
          
          open(5,file=
      1        dir2(1:len2)//'/'//xname6(1:nm6)//'.dat',
      1        status='old')                                 
 c
          read (5,1090) txt
          read (5,1090) txt
          read (5,1125) ifunc,mfunc
          if(ifunc.eq.1)then
            read (5,1090) txt
            read (5,1100) wvlth2,refr1,refi1
            read (5,1090) txt
            read (5,1100) rmin1,rmax1
            read (5,1090) txt
            read (5,1100) delr1,delx1
            read (5,1090) txt
            read (5,1100) rc,cc,xnu       
          elseif(ifunc.eq.2)then
            read (5,1090) txt
            read (5,1100) wvlth2,refr1,refi1
            read (5,1090) txt
            read (5,1100) rmin1,rmax1
            read (5,1090) txt
            read (5,1100) delr1,delx1
            read (5,1090) txt
            read (5,1100) xalf,xgam,xa,xb 
          elseif(ifunc.eq.3)then
            read (5,1090) txt
            read (5,1100) wvlth2,refr1,refi1
            read (5,1090) txt
            read (5,1100) rmin1,rmax1
            read (5,1090) txt
            read (5,1100) delr1,delx1
            read (5,1090) txt
            read (5,1100) rg1,sig1,xnum1
          elseif(ifunc.eq.4)then
            read (5,1090) txt
            read (5,1100) wvlth2,refr1,refi1,refr2,refi2
            read (5,1090) txt
            read (5,1100) rmin1,rmax1,rmin2,rmax2
            read (5,1090) txt
            read (5,1100) delr1,delx1,delr2,delx2
            read (5,1090) txt
            read (5,1100) rg1,sig1,xnum1,rg2,sig2,xnum2
          elseif(ifunc.eq.5)then
            read (5,1090) txt
            read (5,1100) wvlth2,refr1,refi1,refr2,refi2,refr3,refi3
            read (5,1090) txt
            read (5,1100) rmin1,rmax1,rmin2,rmax2,rmin3,rmax3
            read (5,1090) txt
            read (5,1100) delr1,delx1,delr2,delx2,delr3,delx3
            read (5,1090) txt
            read (5,1130) rg1,sig1,xnum1,rg2,sig2,xnum2,rg3,sig3,xnum3 
          endif
          read (5,1090) txt
          read (5,1100) r11,r22,r33,r44,reff,veff
          read (5,1090) txt
          read (5,1095) ccn,bsr,salb,asf,qs,qt
 c        convert the wavelength in cm
          wvlth2=wvlth2*1.0d-4
 c
          dw=abs(wvlth2-wvlth)
 cbaf         if(dw .ge.1.0e-12)then
          if(dw .ge.2.0e-8)then
 c            write(6,1105)wvlth,wvlth2
              write(*,*) "Wavelength mismatch ",dw,wvlth,wvlth2
              stop
          endif
 c
 c....read the prevoiusly computed t matrix
          do i=1,1801
             read(5,1120)(t(i,j),j=1,4),thd(i)
          enddo
 c 
       endif
       if(iref.eq.2)then
       open(18,
      1     access='sequential',form='unformatted',status='scratch')
       open(19,
      1     access='sequential',form='unformatted',status='scratch')
 c
 c      open(20,file=
 c     1     dir4(1:len4)//'/'//'ocn_refl_wl'//cilm//'x'//cwind//'.dat',
 c     2     access='sequential',form='unformatted',status='old')
 c
       open(20,file=
      1     dir4(1:len4)//'/'//'rufrtwl'//cilm//'x'//cwind//'solz00.dat',
      2     access='sequential',form='unformatted',status='old')
 c
 c         read the rough ocean inputs
           rewind(20)
           read (20) xzx,izx,athcnd,aphcnd
           write(*,*)
      1     dir4(1:len4)//'/'//'ocn_refl_wl'//cilm//'x'//cwind//'.dat'
           write(*,901)xzx
           write(*,902)izx
           write(*,903)athcnd
           write(*,904)aphcnd
 901       format('xzx',3f8.3)
 902       format('izx',2i5)
 903       format('athcnd',10f7.1)
 904       format('aphcnd',10f7.1)
           mtha = izx(1)
           mphi = izx(2)
           write (19) (xzx(i),i=1,3),mtha,mphi
 c
           do l=1,mtha
              read (20) thetap,ftrx,ptix
 c            write (18) thetap,(((ftrx(i,j,k),i=1,16),j=1,mtha),
 c    1                            k=1,mphi)
              do k=1,mphi
                 do j=1,mtha
                    do i=1,16
                       txx(i,k,l,j)=ftrx(i,j,k)
                    enddo
                 enddo
              enddo
 c
 c            write (19) (((ptix(i,j,k),i=1,8),j=1,mtha),k=1,mphi)
           enddo
           rewind(18)
           rewind(19)
           rewind(20)
       endif
 c
 c     if iref=3 then read the brdf of the ground surface
 c
       if(iref.eq.3)then
          open(17,file=dir0(1:len0)//'/'//'brdf.dat',
      1        status='old')
 c  
          read(17,'(a)')temp
          read(17,*)nsolbrdf,nthbrdf,nphbrdf,nwbrdf
          do is=1,nsolbrdf
             read(17,160)solbrdf(is)
 160         format(t6,f4.1)
             read(17,165)(wbrdf(i),i=1,nwbrdf)
 165         format(t22,f5.1,t35,f5.1,t48,f5.1,t61,f5.1,
      1             t74,f5.1,t87,f5.1)
 c           check the input wavelength against the wavelengths
 c           in wbrdf            
             mwbrdf=0
             do i=1,nwbrdf
                diffwbrdf=dabs(wvlth*1.0d7-wbrdf(i)) 
                if(diffwbrdf.lt.0.1d-3)mwbrdf=i
             enddo
             if(mwbrdf.eq.0)then
 c              write(6,*)'input wavelength does not matches',
 c    1                   ' with the wavelengths in the buffer',
 c    2                   ' wbrdf'
 c              write(6,*)'input wavelength',wvlth*1.0d7,' nm'
 c              write(6,*)'wavelengths in the buffer wbrdf:',wbrdf
                call closeunits
                stop
             endif
             do j=1,nphbrdf
                m=nphbrdf-j+1
                do k=1,nthbrdf
                   read(17,170)thbrdf(k),phbrdf(j),
      1                        (brdfy(l),l=1,nwbrdf)
 170               format(f6.2,2x,f6.2,10d13.4)
                   brdfx(is,k,m)=brdfy(mwbrdf)
                enddo
             enddo
          enddo
 c
       endif
 c
 c     initialize the buffer & fill the header block 
 c
       if(ipol.eq.0)rho=0.0d0
 c
       do i=1,4625
          bfr1(i)=-777.0
          bfr2(i)=-777.0
          bfr3(i)=-777.0
       enddo
 c
       bfr1(1)=xyear
       bfr1(2)=xmonth
       bfr1(3)=xday
       bfr1(4)=-777.0
       bfr1(5)=wvlth*1.0d4
       bfr1(6)=refr1
       bfr1(7)=refi1
       bfr1(8)=rg1
       bfr1(9)=sig1
       bfr1(10)=rmin1
       bfr1(11)=rmax1
       bfr1(12)=delx1
       bfr1(13)=reff
       bfr1(14)=r11
       bfr1(15)=r22
       bfr1(16)=r33
       bfr1(17)=r44
       bfr1(18)=salb
       bfr1(19)=asf
       bfr1(20)=qs
       bfr1(21)=qt
       bfr1(22)=xzx(3)
       bfr1(23)=xzx(1)
       bfr1(24)=xzx(2)
       bfr1(25)=tr
       bfr1(26)=tm
       bfr1(27)=ta
       bfr1(28)=twat
       bfr1(29)=tautot
       bfr1(30)=rho
       bfr1(31)=ccn
       bfr1(32)=bsr
       bfr1(33)=float(ipol)
       bfr1(34)=float(ipsudo)
       bfr1(41)=dfloat(ifoam)
       bfr1(42)=dfloat(iwatr)
 c     bfr1(43)=chl(iconc)
       bfr1(44)=albwat(ilm)
       bfr1(45)=dfloat(iref)
       bfr1(46)=psrfc
       bfr1(47)=dfloat(ifc)
       bfr1(48)=dfloat(ifunc)
       bfr1(49)=dfloat(mfunc)
       goto(210,220,230,230,230),ifunc
 210   continue
       bfr1(51)=rc
       bfr1(52)=cc
       bfr1(53)=xnu
       goto 260
 220   continue
       bfr1(54)=xalf
       bfr1(55)=xgam
       bfr1(56)=xa
       bfr1(57)=xb
       goto 260
 230   continue
       bfr1(58)=xnum1
       if(ifunc.eq.3)goto 260
       bfr1(59)=rg2
       bfr1(60)=sig2
       bfr1(61)=xnum2
       bfr1(62)=delx2
       bfr1(63)=refr2
       bfr1(64)=refi2
       bfr1(65)=rmin2
       bfr1(66)=rmax2
       if(ifunc.eq.4)goto 260
       bfr1(67)=rg3
       bfr1(68)=sig3
       bfr1(69)=xnum3 
       bfr1(70)=delx3
       bfr1(71)=refr3
       bfr1(72)=refi3
       bfr1(73)=rmin3
       bfr1(74)=rmax3
 260   continue
       bfr1(76)=itrans  !flag for computing diffused transmiss
 c
 c*****format statements******************************************************
  100  format(2i4,1x,1p6e11.3)                                           
  125  format(a)
  126  format(6d12.4,i4)
  128  format(4x,6d12.4)
  145  format(i5)
  150  format(5x,1p6e12.4)                                               
   195 format (1x,'xzx',1p10d11.2)
   196 format (1x,'izx',2i5)
   197 format (1x,'athcnd',1p10d11.2)
   198 format (1x,'aphcnd',1p10d11.2)
   199 format (1x,'thetap',1p10d11.2)
   200 format (1x,'ftrx',1p10d11.2)
 c  205 format('xzx1-3,mtha,mphi'/1p3e11.3,0p2i4)
   201 format (1x,'pti',1p10d11.2)
   505 format (1x,'end of file encountered , check the program')
  1090 format(a)
  1095 format(1p4e11.3,1p2e12.4)
  1100 format(1p7e11.3)
  1105     format(1x,'input wavelength and phase matrix wavelength'
      1       ' do not match.'/1x,'input wavelength=',1pe12.4,
      2         'phase matrix wavelength=',1pe12.4)
  1120 format(1p2d18.12,1p2d19.12,0pf6.2)
  1125 format(i5,11x,i5)
  1130 format(9f8.5)
  1160 format(a)
 c****************************************************************************
       return                                                            
       end                                                               
 c***********************************************************************

index

an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now

closeunits

subroutine closeunits

Definition: closeunits.f:2

ccn

subroutine ccn(ifunc, irgm)

Definition: phs.f:940

readin

subroutine readin(ilm, irh, isd, itau)

Definition: readin.f:2

real

#define real

Definition: DbAlgOcean.cpp:26

float rc(float x, float y)

abs

#define abs(a)

Definition: misc.h:90

convtc

subroutine convtc(num, nchar, loc)

Definition: ocn.f:156