fltocn_new.f

Go to the documentation of this file.

      subroutine fltocn_new
c     subroutine fltocn(fio,bmu,dmu,dmus2,eo,ddphi,amuo,pi,sumc,sumcpi,
c    1                  sumdwn,calb,rmu,conv,eox,kkx,jjjj,nmum1,jpart)
c
c     treat the lower boundary of the atmosphere as flat ocean and 
c     reflect all incident light (diffuse and direct) according to 
c     fresnel law
c
c***********************************************************************
c.....include the common and declaration statemnets.....................
      implicit real*8 (a-h,o-z)
c
      include 'common_all.cmn'
      real*8 zx(4,4),temp(4,50,46)
c***********************************************************************
 
c***********************************************************************
c     implicit real*8 (a-h,o-z)
c     real*8 fio(4,50,46),bmu(51),dmu(51),dmus2(51),eo(4)
c     real*8 temp(4,50,46),z(4,4),tto(4),ttf(4),tth(4),ttd(4)
c     real*8 zx(4,4),rmu(51),eox(4)
c***********************************************************************
c
c     write(*,*)'welcome to fltocn'
      xzx=-1.0d-8                                                       
      xr=1.3340d0
      xi=0.0d0
c     initalize the glint buffer
      do i=1,nmum1
         do  j=1,jpart
             do  k=1,4
                 fglint(k,i,j)=0.0d0
             enddo 
         enddo
      enddo
      do 500 it=1,jjjj
         itp=(nmum1-it)+1
c      write(*,*)'ready to begin the 400 loop'
         cost=bmu(it)
         thetac= dacos(cost)                                               
         sinsq = (1.0d0-cost**2)                                           
         sinsq2 = sinsq**2                                                 
         a=(xr**2-xi**2-1.0d0+cost**2)                            
         b=(-2.0d0*xr*xi)                                        
         r=dsqrt(a**2+b**2)                                        
         tmr=cost**2*r                                                
         xp=dsqrt(2.0d0*r+2.0d0*a)                                  
         xmm=(2.0d0*r-2.0d0*a)                                  
         if(xmm .lt. 0.0d0 .and. xmm .gt. xzx) xmm=-xmm                    
         xm=dsqrt(xmm)                                                     
         dnmr1=(sinsq2+tmr+cost*sinsq*xp)                           
         qrmu=(sinsq2-tmr)/dnmr1                                     
         qimu=(cost*sinsq*xm)/dnmr1                                    
         dnmr2=(cost**2+r+cost*xp)                                  
         rrr=(cost**2-r)/dnmr2                                      
         rri=(cost*xm)/dnmr2                                          
         rer=qrmu*rrr-qimu*rri                                      
         rei=qimu*rrr+qrmu*rri                                      
         r11=rer**2+rei**2                                          
         r22=rrr**2+rri**2                                          
         r33=rer*rrr+rei*rri                                        
         r34=rri*rer-rei*rrr                                        
         zx(1,1)=r11
         zx(2,2)=r22
         zx(3,3)=r33
         zx(3,4)=r34
         zx(4,3)=-r34
         zx(4,4)=r33
         
c        write(6,120)cost,r11,r22,r33,r34
120      format(f8.4,1p4e15.5)
c        
c
c        write(*,*)'zx matrix',zx
c        apply the zx matrix to the incident radiation
c        solid angle of the incident beam=solid angle of the reflected beam
         do ip=1,jpart
            do i=1,4
               sum2=0.0d0
               do j=1,4
                  sum2=sum2+zx(i,j)*fio(j,it,ip)
               enddo
               temp(i,itp,ip)=sum2
            enddo
         enddo
         if(it.eq.kkx)then
c           write(*,*)'eo',eo
            sang=dmu(kkx)*ddphi
            do i=1,4
               sum1=0.0d0
               do j=1,4
                  sum1=sum1+zx(i,j)*eo(j)/sang
               enddo
               temp(i,itp,1)=temp(i,itp,1)+sum1
               fglint(i,itp,1)=sum1
            enddo
         endif
 500  continue
c
c     tranfer the stokes parameter from the temp buffer to fio
c     only for the upward direction
      do it=jjjj+1,nmum1
         do ip=1,jpart
            do i=1,4
               fio(i,it,ip)=temp(i,it,ip)
            enddo
         enddo
      enddo
c
c     compute the albedo of the flat ocen
      sumc=0.0d0
      do it=jjjj+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
c     write(*,*)'sumc,sumcpi,sumdwn,sumd',sumc,sumcpi,sumdwn,sumd
      calb=sumc/sumd
c     write(*,*)'ocean albedo=',calb 
      return                                                            
      end                                                               
c***********************************************************************