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 'common_all.cmn'
c
c**************************************************************************
c
c     compute the downward diffuse radiation at each level
c
      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,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))+
     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,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))+
     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=nolyrp)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=nolyrp)fio
      endif
c
      if(iref.eq.3)then
         call brdfg(fio,brdfx,bmu,dmu,dmus2,eo,ddphi,amuo,pi,sumc,
     1               sumcpi,sumdwn,calb,kkx,jjjj,nmum1,jpart)
        write(iwrt,rec=nolyrp)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(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))+
     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(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))+
     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************************************************************************