rsea_new.f

Go to the documentation of this file.

      subroutine rsea_new
c
c     treat the lower boundary of the atmosphere as non lambertian 
c     surface (rough ocean surface) and reflect all incident light
c     (diffuse and direct light) 
c     using the reflection matrix of the surface
c
c***********************************************************************
      implicit real*8 (a-h,o-z)
      include 'common_all.cmn'
      real*8 fiit(4,50,46),tempa(4,50,46),tempb(4,50,46),tempc(4,50,46)
      real*8 tss(4,50,46),temph(4,50,46),tempg(4,50,46)
      real*8 glint_diff(4,50,46),tempd(4,50,46)
c***********************************************************************
c
c     write(*,*)'welcome to rsea_new'
c     initialize the buffers
      do i=1,nmum1
         do j=1,jpart
            do k=1,4
               tempa(k,i,j)=0.0d0
               tempb(k,i,j)=0.0d0
               tempc(k,i,j)=0.0d0
               tempd(k,i,j)=0.0d0
               tempg(k,i,j)=0.0d0
               temph(k,i,j)=0.0d0
               tss(k,i,j)=0.0d0
            enddo
         enddo
      enddo
c     total downward flux
      call fluxlvl(fio,sumdnz,0)
      pflx0=(eo(1)+eo(2))*bmu(kkx)
      pflx=pflx0+sumdnz
c
c     write(*,*)'pflx0,sumdnz,pflx',pflx0,sumdnz,pflx
 
c
c     contribution from the direct beam
      do it=1,jjjj
         itt=nmum1-it+1
         do ip=1,jpart
            do is=1,4
               do k=1,4
                  ms=(is-1)*4+k
                 tempa(is,itt,ip)=tempa(is,itt,ip)+
     1                            txx(ms,ip,kkx,it)*eo(k)
               enddo
            enddo
         enddo
      enddo
c
c
      call fluxlvl(tempa,flxtempa,1)
c           print the upwelling radiance due to dir. flux
c           do i=jjj,nmum1
c              call radnce(pi,conv,bmu,the,tempa,i,jpart,jphi)
c           enddo
c
c     contribution from diffuse beam (integrate over thetaprime & phiprime)
      do it=1,jjjj
         itt=nmum1-it+1
         do ip=1,jpart
          do is=1,4
            do itp=1,jjjj
               sangi=dabs(dmu(itp))*ddphi
               do ipp=1,nophi
                  jpp=ipp
                  jpx=iabs(ipp-ip)+1
                  if(ipp.gt.jpart)jpp=nophi-jpp+2
                  if(jpx.gt.jpart)jpx=nophi-jpx+2
                     do i=1,2
                        fiit(i,itp,jpp)=fio(i,itp,jpp)
                        fiit(i+2,itp,jpp)=fio(i+2,itp,jpp)
                        if(ipp.gt.jpart)fiit(i+2,itp,jpp)=
     1                                  -fiit(i+2,itp,jpp)
                     enddo
                  prodt3a=0.0d0
                  do k=1,4
                     m=(is-1)*4+k
                     prodt3a=prodt3a+txx(m,jpx,itp,it)*fiit(k,itp,jpp)
                  enddo
                  tempc(is,itt,ip)=tempc(is,itt,ip)+prodt3a*sangi
               enddo
            enddo
          enddo
         enddo
      enddo  
c
 
c     save the glint contribution from the direct beam 
c     into fglint buffer
      do i=1,nmum1
         do j=1,jpart
            do k=1,4
               fglint(k,i,j)=tempa(k,i,j)
            enddo
         enddo
      enddo
c
      call fluxlvl(tempc,flxtempc,1)
c     write(*,*)'flux_tempc=',flxtempc
c
      do i=jjj,nmum1
         ii=i-jjjj
         sang=dabs(dmus2(i))*ddphi
         do j=1,jpart
            do k=1,4
               tss(k,i,j)=tempa(k,i,j)+tempc(k,i,j)
               tss(k,ii,j)=0.0d0
            enddo
         enddo
      enddo
c
c     compute the upwelling fluxes due to each component
      call fluxlvl(tempa,flxza,1)
      call fluxlvl(tempc,flxzc,1)
      albdtot=(flxza+flxzc)/pflx
      albrdr(ksza)=flxza/pflx
      albrdf(ksza)=flxzc/pflx
c
c     correct the upwelling radiance for foam and underwater radiances
c     foam correction
c
      call foam
c
c     underwater contribution
c     compute the upwelling radiance from the albedo of water just
c     below the ocean surface
      if(iwatr.eq.1)then
        call watlev
      endif
c     compute the flux of the water leaving radiance
      qsumzz=0.0d0
      do i=jjj,nmum1
         qsumz=0.0d0
         do j=1,jpart
            if(j.eq.1 .or. j.eq.jpart)then
               qsumz=qsumz+radxi(i,j)
            else
               qsumz=qsumz+2.0d0*radxi(i,j)
            endif
         enddo
         qsumzz=qsumzz+qsumz*dabs(dmus2(i))*ddphi
      enddo
      albwl(ksza)=qsumzz
c     write(*,176)qsumzz
176   format('flux of the waterleavin radiance',1pe12.4)
c
      do it=jjj,nmum1
        do ip=1,jpart
         do ks=1,4
           if(ks.le.2)then
             fio(ks,it,ip)=(1.0d0-fracfm)*tss(ks,it,ip)+
     1                    0.5d0*xifm+0.5d0*radxi(it,ip)
           else
             fio(ks,it,ip)=(1.0d0-fracfm)*tss(ks,it,ip)
           endif
         enddo
        enddo
      enddo
c
      do it=jjj,nmum1
         itp=nmum1-it+1
         do ip=1,jpart
               raddir(ksza,itp,ip)=tempa(1,it,ip)+tempa(2,it,ip)
               radocn(ksza,itp,ip)=radxi(it,ip)
               radsky(ksza,itp,ip)=tempc(1,it,ip)+tempc(2,it,ip)
         enddo
      enddo
c
      call fluxlvl(fio,sumc,1)
c     write(*,*)'flux above the ocean surface',sumc
c     compute the albedo of the surface
      sumcpi=sumc*pi
      sumd=pflx
c     write(*,*)'sumc,sumcpi,sumdwn,sumd',sumc,sumcpi,sumdwn,sumd
      calb=sumc/sumd
c     write(*,*)'surface albedo=',calb 
      return                                                            
      end                                                               
c***********************************************************************