main.f

Go to the documentation of this file.

c***********************************************************************
c                                                                       
c                                                                       
c     radiative transfer program (part 2)                               
c                                                                       
c                                                                       
c***********************************************************************
c.....include the common and declaration statemnets.....................
      implicit real*8 (a-h,o-z)
c
      real*8 glint_tmp(4,50,46),nrhum,nsd
      include 'common_all.cmn'                                         
c***********************************************************************
c                                                                       
c
c     get the upper and lower indices for wavelength, size dist. and 
c     optical thickness.  also, flags for wind speed, foam, water leaving
c     radiances and chlorophyll conc.
c
      nrhum = 8
      nsd = 80
      izia=1
c      iglint=1         !subtract the glint radiance
      iocn=1           !wind speed in sigma of slope dist.
      call inusr
      if(iref.lt.2)then
        iwnd1=1
        iwnd2=1
      endif
c
      do 878 ilm=ilm1,ilm2,1
      do 877 isd=isd1,isd2,1
c      do 874 irh=irh1,irh2,1
      do 876 itau=itau1,itau2,1
      do 875 iwnd=iwnd1,iwnd2,1
 
      irh = (isd-1)*nrhum/nsd + 1
      write(*,*) ilm, isd, irh, itau, iwnd
c
c     open files and read input created in part 1 of the program        
      iwind=iwnd
      call readin(ilm,irh,isd,itau)
c
      if(icrft.eq.1)then
         call crftlvl
      endif
c
      xlam=wvlth                                                        
      alw=albwat(ilm)
c
c     loop over solar zenith angle                                      
      nsza=0                                                            
      do 5566 ksza=ithe01,ithe02,1
      isza=ksza-ithe01+1
      msza(isza)=ksza
      nsza=nsza+1                                                       
      call iniclz                                                       
      call angl                                                         
      if(isza.eq.1)then
         call hdrmds(ithe01,ithe02)
      endif
c
c     if ipsudo=1 then determine optical path length for                
c     spherical atmosphere                                              
      if(ipsudo.eq.1)then                                               
         call spathz(amuo,htp,taur,taum,tauabs,totsp,nmodl)             
      endif                                                             
c
c     determine an average value of phase matrix in the forward/        
c     backward direction                                                
      if(ifc.eq.1)then                                
             call hump (const,t,pp,nmum1,lphi,rmu,thd,0,nmu)  
             call hump (const,t,qq,nmum1,lphi,rmu,thd,1,nmu) 
      endif                                                             
c
c     print summry of the input parameters                              
c     call summry                                                       
c
c     if iref=2 then read the rough ocean input created in              
c     part 1 of the program                                             
c
c     if(itrans.eq.1)then
c        call anglw
c     endif
c
      if(iref .eq.2)then                                                
         read(19)xrw,xiw,vz,mtha,mphi 
         call anglw
      endif                                                             
c                                                                       
c     compute the attenuation parameters                                
      call attenew
c                                                                       
c     compute the scattering phase matrices (rayleigh/aerosol)          
      call phase                                                        
c
c     begin radiative transfer calculations                             
c
      ibgn=1
      iend=1
      if(nsza.eq.1 .and. (iref.eq.0 .or. itrans.eq.1) )then
        iend=2
      endif
c
c     illumination from the top of the atmosphere(kzz=1)                
c     illumination from the bottom of the atmosphere(kzz=2)             
c
      do 20 kzz=ibgn,iend,1
            jpass=1
            nump=npass1
            call iniclz2
            if(nsza.eq.1 .and. (iref.eq.0 .or. itrans.eq.1) )then
               nump=npass2
            endif
            if(kzz.eq.1)then
               call single
               minitr=4
            endif
            if(kzz.eq.2 .and. iref.eq.0)then
               call grefl 
               call snglup
               minitr=4
            endif
            if(kzz.eq.2 .and. itrans.eq.1)then
               call focn_below
               call snglup
               minitr=4
            endif
c
            call getrad
c              
         do 15 ka=1,nump                                           
c           initialize the disk units
            kdx=ka-(ka/2)*2
            if(kdx.eq.1)then
              irad=64
              iwrt=65
            else
              irad=65
              iwrt=64
            endif
            if(kzz.eq.1)then
              call multp1d
            elseif(kzz.eq.2 .and. iref.eq.0)then
              call multp2d
            elseif(kzz.eq.2 .and. itrans.eq.1)then
              call multp2d_trans
            endif
            jpass=jpass+1
            call getrad
            call contst
            call fiotb
            if(iactflx.eq.1)then
               call actflx
            endif
            if(jpass.ge.minitr)then
               if(d3.le.0.1 .or. jpass.ge.20)goto 17
            endif
 15         continue                                                    
 17         continue
c
            if(kzz.eq.1 .and. iglint.eq.1)then
              do k=1,4
                 do i=jjj,nmum1
                    do j=1,jpart
                       glint_tmp(k,i,j)=fglint(k,i,j)*
     1                                  dexp(-tautot/dabs(bmu(i)))
                    enddo
                 enddo
              enddo 
c             if(kzz .eq.2)then
c             do i=jjj,nmum1
c                 call radnce(pi,conv,bmu,the,fio,i,jpart,jphi)
c             enddo
c             endif
            endif
c                                                                       
            call geocor
            if(icrft.eq.1)then
               call crftgcr
            endif
c
            if(iactflx.eq.1)then
               call actfgrc
            endif
c
            if(isurf.eq.1)then
               call surfgcr
            endif
c
c           illumination brom botom is done only once
            if(kzz.eq.2)goto 30           
20    continue                                                          
c
30    continue
c
      rewind 21                                                         
      if (ifc.eq.1) then
         rewind 5
      endif
      rewind 11                                                         
      rewind 12                                                         
      rewind 18                                                         
      rewind 19                                                         
c
5566  continue                                                          
c
c     create output datasets for upwelling and downwelling radiances    
c     leaving the top and bottom of the atmosphere                      
c           
      if(iref.eq.1 .and. itrans.eq.1)then
        call outdt_trans
      endif
c                                                            
      if((iref.eq.1 .or. iref.eq.2 .or. iref.eq.3) .and.itrans.eq.0)then
        call outdty
        if(icrft.eq.1)then
           call outcrfty
        endif
        if(iactflx.eq.1)then
           call outactfy
        endif
        if(isurf.eq.1)then
           call outsurfy
           if(iref.eq.2)then
c             call outsurfdir
              call outsurfocn
c             call outsurfsky
           endif
        endif
      endif
c
      if(iref.eq.0 .and.itrans.eq.0)then
        call outdtz
        if(icrft.eq.1)then
           write(*,*)'ready to call outcrftz'
           call outcrftz
        endif
        if(iactflx.eq.1)then
           call outactfz
        endif
      endif
c                                                                       
c     close the temp files
c
      call closeunits
      call closemore
c
875   continue
876   continue                                                          
874   continue
877   continue                                                          
878   continue                                                          
c                                                                       
c********************************************************************** 
c                   format statements                                   
355    format( t10,'fresnel reflection (by a rough surface)',1x,
     1 'at the lower boundary '/ t10,'refractive index',t45,'=',1pe15.5,
     2 '-',1pe15.5,'i' / t10,'velocity',t45,'=', 1pe15.5,'meter/sec'/)
c***********************************************************************
      end                                                               
c**********************************************************************