brdfg.f

Go to the documentation of this file.

      subroutine brdfg(fio,brdfx,bmu,dmu,dmus2,eo,ddphi,amuo,pi,sumc,
     1                  sumcpi,sumdwn,calb,kkx,nx,nmum1,jpart)
c
c     treat the lower boundary of the atmosphere as non lambertian 
c     surface and reflect all incident light (diffuse and direct) 
c     using the brdf of the surface
c
c***********************************************************************
      implicit real*8 (a-h,o-z)
      real*8 fio(4,50,46),cosmu(51),dmu(51),dmus2(51),eo(4)
      real*8 temp(4,50,46),z(4,4),tto(4),ttf(4),tth(4),ttd(4)
      real*8 brdfx(25,25,46)
c***********************************************************************
c
      write(*,*)'welcome to brdfg'
      do i=1,(nx-1)
         do j=1,jpart
c           write(*,100)i,j,(fio(k,i,j),k=1,4)
100         format('fio',i3,i3,1p4e12.3)
         enddo
      enddo
c
c     write(*,*)'cosmu',cosmu
c     write(*,*)'dmu',dmu
c     write(*,*)'dmus2',dmus2
c     write(*,*)'eo',eo
c     write(*,*)'ddphi,amuo',ddphi,amuo
c     write(*,*)'(nx-1),nmum1,jpart',(nx-1),nmum1,jpart
c
      mph2=2*(jpart-1)
      do it=1,(nx-1)
         itt=nmum1-it+1
         do ip=1,jpart
c                             
            sumtp=0.0d0
            do itp=1,(nx-1)
               sangi=dmu(itp)*ddphi
               do ipp=1,mph2
                  ippz=ipp
                  if(ippz.gt.jpart)ippz=mph2-ipp+2
                  ipx=iabs(ipp-ip)+1
                  if(ipx.gt.jpart)ipx=mph2-ipx+2
                  fdif=(fio(1,itp,ippz)+fio(2,itp,ippz))*
     1                  sangi*cosmu(itp)
                  sumtp=sumtp+brdfx(itp,it,ipx)*(1/pi)*fdif
                  if(itp.eq.kkx .and. ipp.eq.1)then
                     fsun=(eo(1)+eo(2))*amuo
                     sumtp=sumtp+brdfx(itp,it,ipx)*(1/pi)*fsun
                  endif
               enddo
            enddo
            do i=1,2
               fio(i,itt,ip)=0.5d0*sumtp
               fio(i+2,itt,ip)=0.0d0
            enddo
         enddo
      enddo
c     compute the albedo of the surface
      sumc=0.0d0
      do it=(nx-1)+1,nmum1
         sum3=fio(1,it,1)+fio(2,it,1)+fio(1,it,jpart)+fio(2,it,jpart)
         do ip=2,jpart-1
            sum3=sum3+2.0d0*(fio(1,it,ip)+fio(2,it,ip))
         enddo
         sumc=sumc+sum3*dabs(dmus2(it))*ddphi
      enddo
      sumcpi=sumc*pi
      sumd=sumdwn+(eo(1)+eo(2))*amuo
      write(*,*)'sumc,sumcpi,sumdwn,sumd',sumc,sumcpi,sumdwn,sumd
      calb=sumc/sumd
      write(*,*)'surface albedo=',calb 
      return                                                            
      end                                                               
c***********************************************************************