multp1d.f

Go to the documentation of this file.

      subroutine multp1d
c
c     compute the multiple scattering contribution when the atmosphere is
c     illuminated from the top.  also, for each level, store the 
c     upwelling/downwelling diffuse radiation.
c
c*************************************************************************
c.....includes the common blocks
      implicit real*8 (a-h,o-z)
      include 'afrt_rt2.cmn'
c
c**************************************************************************
c
c     compute the downward diffuse radiation at each level
c
      do i=1,4
      do j=1,2*nsz
      do k=1,nph
      sngla(i,j,k) = 0.0
      snglb(i,j,k) = 0.0
      enddo
      enddo
      enddo
 
      isng=53
      read(irad,rec=1)fiib
      write(iwrt,rec=1)fiib
      do il=1,nolyr
         ilp=il+1
         ilm=il-1
          if(il.eq.1)then
             read(iwrt,rec=il)fiib
             read(irad,rec=ilp)fio
             read(isng,rec=il)sngla
             read(isng,rec=ilp)snglb
c        determine average intensity at the center of the layer il
             do i=1,nmum1
             do j=1,jpart
                do k=1,4
                   fiic(k,i,j)=0.50d0*(fio(k,i,j)+fiib(k,i,j))
                   ftmp(k,i,j)=0.0d0
                enddo
             enddo
             enddo
             tmsl=dtmm(il)*qsqt*const
             trsl=dtrr(il)*conr
             dlyr=dtot(il)
             if(ifc.eq.0)tmsl=0.0d0
             call mdiffn(1,nx-1,il,dlyr)
             do i=1,nx-1
                do j=1,jpart
                   do k=1,4
                      fio(k,i,j)=fiib(k,i,j)*emdtm(il,i)+
     1                       ftmp(k,i,j)*(1.0d0-emdtm(il,i))+
     2                       snglb(k,i,j)-sngla(k,i,j)*emdtm(il,i)
                   enddo
                enddo
             enddo
             write(iwrt,rec=ilp)fio
          else
             read(iwrt,rec=ilm)fiib
             read(isng,rec=ilm)sngla
             read(isng,rec=ilp)snglb
             do i=1,nmum1
                 do j=1,jpart
                    do k=1,4
                       fiic(k,i,j)=fio(k,i,j)
                       ftmp(k,i,j)=0.0d0
                    enddo
                 enddo
             enddo
             read(irad,rec=ilp)fio
             tmsl=(dtmm(ilm)+dtmm(il))*qsqt*const
             trsl=(dtrr(ilm)+dtrr(il))*conr
             dlyr=(dtot(ilm)+dtot(il))
             if(ifc.eq.0)tmsl=0.0d0
             call mdiffn(1,nx-1,il,dlyr)
             do i=1,nx-1
                do j=1,jpart
                   do k=1,4
                      fio(k,i,j)=fiib(k,i,j)*emdtm(ilm,i)*emdtm(il,i)+
     1              ftmp(k,i,j)*(1.0d0-emdtm(ilm,i)*emdtm(il,i))+
     2              snglb(k,i,j)-sngla(k,i,j)*emdtm(ilm,i)*emdtm(il,i)
                   enddo
                enddo
             enddo
             write(iwrt,rec=ilp)fio
c
          endif
      enddo
c
      if(iref.eq.1)then
         call fltocn_new
         write(iwrt,rec=nolyr+1)fio
      endif
 
      if(iref.eq.2)then
c       call rsea_new
        if(iocn.eq.1) then
           call rsea_new
        else
c          call RSEA
        endif
        write(iwrt,rec=nolyr+1)fio
      endif
c
      if(iref.eq.3)then
         call brdfg(fio,brdfx,cosmu,dcmu,dcmusq,eo,ddphi,amuo,pi,sumc,
     1               sumcpi,sumdwn,calb,kkx,nx-1,nmum1,jpart)
        write(iwrt,rec=nolyr+1)fio
      endif
c
c     compute the upward diffuse radiation at each level
c
      do il=1,nolyr
         im=nolyr-il+1
         imp=im+1
         impp=imp+1
         if(il.eq.1)then
             read(iwrt,rec=imp)fiib
             read(iwrt,rec=im)fio
             read(isng,rec=imp)sngla
             read(isng,rec=im)snglb
             do i=1,nmum1
                 do j=1,jpart
                    do k=1,4
                       fiic(k,i,j)=0.50d0*(fio(k,i,j)+fiib(k,i,j))
                       ftmp(k,i,j)=0.0d0
                    enddo
                 enddo
             enddo
             tmsl=dtmm(im)*qsqt*const
             trsl=dtrr(im)*conr
             dlyr=dtot(im)
             if(ifc.eq.0)tmsl=0.0d0
             call mdiffn(nx,nmum1,im,dlyr)
             do i=nx,nmum1
                do j=1,jpart
                   do k=1,4
                      fio(k,i,j)=fiib(k,i,j)*emdtm(im,i)+
     1                       ftmp(k,i,j)*(1.0d0-emdtm(im,i))+
     2                       snglb(k,i,j)-sngla(k,i,j)*emdtm(im,i)
                   enddo
                enddo
             enddo
             write(iwrt,rec=im)fio
         else
             read(iwrt,rec=impp)fiib
             read(isng,rec=impp)sngla
             read(isng,rec=im)snglb
             do i=1,nmum1
                 do j=1,jpart
                    do k=1,4
                       fiic(k,i,j)=fio(k,i,j)
                       ftmp(k,i,j)=0.0d0
                    enddo
                 enddo
             enddo
             read(iwrt,rec=im)fio
             tmsl=(dtmm(imp)+dtmm(im))*qsqt*const
             trsl=(dtrr(imp)+dtrr(im))*conr
             dlyr=(dtot(imp)+dtot(im))
             if(ifc.eq.0)tmsl=0.0d0
             call mdiffn(nx,nmum1,im,dlyr)
             do i=nx,nmum1
                do j=1,jpart
                   do k=1,4
                      fio(k,i,j)=fiib(k,i,j)*emdtm(imp,i)*emdtm(im,i)+
     1                  ftmp(k,i,j)*(1.0d0-emdtm(imp,i)*emdtm(im,i))+
     2                  snglb(k,i,j)-sngla(k,i,j)*emdtm(imp,i)*emdtm(im,i)
                   enddo
                enddo
             enddo
             write(iwrt,rec=im)fio
c
         endif
      enddo
c
      return
      end
c
c************************************************************************