multp2d.f

Go to the documentation of this file.

      subroutine multp2d
c
c     compute the multiple scattering contribution when the atmosphere is
c     illuminated from the bottom.  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 'common_all.cmn'
c
c**************************************************************************
c
c     compute the downward diffuse radiation at each level
c
      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
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,jjjj,il,dlyr)
         do i=1,jjjj
            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))
               enddo
            enddo
         enddo
         write(iwrt,rec=ilp)fio
      else
         read(iwrt,rec=ilm)fiib
         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,jjjj,il,dlyr)
         do i=1,jjjj
            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))
               enddo
            enddo
         enddo
         write(iwrt,rec=ilp)fio
      endif
      enddo
c
c     compute the upward diffuse radiation at each level
      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
         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(jjj,nmum1,im,dlyr)
         do i=jjj,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))
               enddo
            enddo
         enddo
         write(iwrt,rec=im)fio
      else
         read(iwrt,rec=impp)fiib
         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(jjj,nmum1,im,dlyr)
         do i=jjj,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))
               enddo
            enddo
         enddo
         write(iwrt,rec=im)fio
      endif
      enddo
c
      return
      end
c***********************************************************************