bri_tem_nonsnow.f90

来自「CLM集合卡曼滤波数据同化算法」· F90 代码 · 共 2,381 行 · 第 1/5 页

! 2nd-order shadowing correctionc
    sheff=shad_cor(i,j)

	xx=spectrum(kx+u,ky+v,fm)
	yy=spectrum(ksx+u,ksy+v,fn)
    
	if(dabs(xx) .le. 0d0 .and. dabs(yy) .le. 0d0)  go to 13
	 
    sp1_LOG=dlog(xx)+dlog(yy)+dlog(w1(i))+dlog(w2(j))+dlog(z1(i))


   
!    sp1=spectrum(kx+u,ky+v,fm)*spectrum(ksx+u,ksy+v,fn) &
!               *w1(i)*w2(j)*z1(i)
!   	if (dabs(sp1).gt.0.d0) then
!	write(*,*) "SD"
!	endif

   sum3hh=sum3hh+funhh(u,v)*sheff*dexp(tempc1_LOG(n)+tempc2_LOG(m)+sp1_LOG)
   sum3vv=sum3vv+funvv(u,v)*sheff*dexp(tempc1_LOG(n)+tempc2_LOG(m)+sp1_LOG)
   sum3hv=sum3hv+funhv(u,v)*sheff*dexp(tempc1_LOG(n)+tempc2_LOG(m)+sp1_LOG)
   sum3vh=sum3vh+funvh(u,v)*sheff*dexp(tempc1_LOG(n)+tempc2_LOG(m)+sp1_LOG)



13        continue
!           if(cdabs(sum3vh-old).le.torlant) go to 16
!            old=sum3vh
14        continue
15     continue
16     continue
!----------------------------------------c
! sum up each term                       c
!----------------------------------------c
      crterm(1)=dreal(dconjg(fhh)*(funhh(-kx,-ky)*sum1 &
                          +funhhs(-ksx,-ksy)*sum2			 &
                          +sum3hh/pi))
      crterm(2)=dreal(dconjg(fvv)*(funvv(-kx,-ky)*sum1 &
                          +funvvs(-ksx,-ksy)*sum2			 &
                          +sum3vv/pi)) 
      crterm(3)=dreal(dconjg(fhv)*(funhv(-kx,-ky)*sum1 &
                          +funhvs(-ksx,-ksy)*sum2			 &
                          +sum3hv/pi))
      crterm(4)=dreal(dconjg(fvh)*(funvh(-kx,-ky)*sum1 &
                          +funvhs(-ksx,-ksy)*sum2			 &
                          +sum3vh/pi))    
17    continue
!--------------------------------------c
! end of computation of cross terms    c
!--------------------------------------c
      if(icort.eq.0) go to 35
!--------------------------------------c
! evaluate  complementary term         c
!--------------------------------------c
      tempr1(1)=(sig*ksz)**2
      tempr2(1)=(sig2*kz*ksz)
      tempr3(1)=(sig*kz)**2
	  tempr1_LOG(1)=dlog(tempr1(1))
	  tempr2_LOG(1)=dlog(tempr2(1))
	  tempr3_LOG(1)=dlog(tempr3(1))
      do 18 n=2,mm
        fn=float(n)
	    tmp_LOG=dlog(fn)
		tempr2_LOG(n)=tempr2_LOG(1)+tempr2_LOG(n-1)-tmp_LOG
		tempr3_LOG(n)=tempr3_LOG(1)+tempr3_LOG(n-1)-tmp_LOG
		tempr1_LOG(n)=tempr1_LOG(1)+tempr1_LOG(n-1)-tmp_LOG
       !tempr1(n)=tempr1(1)*tempr1(n-1)/fn
18     continue
!       do 19 n=2,mm
!       fn=float(n)
       !tempr2(n)=tempr2(1)*tempr2(n-1)/fn
19     continue
!do 20 n=2,mm
!      fn=float(n)
!       tempr3(n)=tempr3(1)*tempr3(n-1)/fn
!20     continue
sum1=0.0
sum2=0.0
sum3=0.0
sum1_LOG=0.0
sum2_LOG=0.0
sum3_LOG=0.0
sum4hh=(0.0,0.0)
sum4vv=(0.0,0.0)
sum4hv=(0.0,0.0)
sum4vh=(0.0,0.0)
sum5hh=(0.0,0.0)
sum5vv=(0.0,0.0)
sum5hv=(0.0,0.0)
sum5vh=(0.0,0.0)
old=(0.0,0.0)
rcom_LOG=-sig2*(ksz*ksz+kz*kz)+dlog(0.125*k*k)
do 21 i=1,mm
   fi=float(i)
   spe=spectrum(ksx-kx,ksy-ky,fi)
!	sum1=sum1+tempr1(i)*spe
!	sum2=sum2+tempr2(i)*spe
!	sum3=sum3+tempr3(i)*spe
!   tmp_LOG=dlog(spe)
sum1_LOG=sum1_LOG+dexp(tempr1_LOG(i)+rcom_LOG)*spe
sum2_LOG=sum2_LOG+dexp(tempr2_LOG(i)+rcom_LOG)*spe
sum3_LOG=sum3_LOG+dexp(tempr3_LOG(i)+rcom_LOG)*spe

21     continue
sum1=sum1_LOG
sum2=sum2_LOG
sum3=sum3_LOG
       if(muti.eq.0) go to 30 
!------------------------------------------c
! calculate the mutiple scattering parts:  c
! two double integrations                  c
!------------------------------------------c      
      do 25 n=1,mm
fn=float(n)
       do 24 m=1,mm
  fm=float(m)
shh2=(0.0,0.0)
svv2=(0.0,0.0)
shv2=(0.0,0.0)
svh2=(0.0,0.0)
shh3=(0.0,0.0)
svv3=(0.0,0.0)
shv3=(0.0,0.0)
svh3=(0.0,0.0)
!---------------------------------c    
do 23 i=1,nw
   do 22 j=1,nw
   u=z1(i)*dcos(z2(j))
   v=z1(i)*dsin(z2(j))
! 2nd-order incident angle
     ql=dasind(dsqrt(u**2+v**2)/k)
!             call shadow(back,ql,ql,sheff)
      sheff=shad_cor(i,j)
! 2nd-order shadowing correctionc 
   spe=spectrum(ksx+u,ksy+v,fn)*spectrum(kx+u,ky+v,fm)
   spcom=spe*w1(i)*w2(j)*z1(i)
     fphh=funhh(u,v)
     fpvv=funvv(u,v)
     fphv=funhv(u,v)
     fpvh=funvh(u,v)
   shh2=shh2+cdabs(fphh)**2*spcom*sheff
   svv2=svv2+cdabs(fpvv)**2*spcom*sheff
   shv2=shv2+cdabs(fphv)**2*spcom*sheff
   svh2=svh2+cdabs(fpvh)**2*spcom*sheff
     km=u+kx+ksx
     kn=v+ky+ksy 
    if((km*km+kn*kn).lt.(k*k)) then
      shh3=shh3+fphh*dconjg(funhh(-km,-kn))*spcom*sheff
      svv3=svv3+fpvv*dconjg(funvv(-km,-kn))*spcom*sheff
      shv3=shv3+fphv*dconjg(funhv(-km,-kn))*spcom*sheff
      svh3=svh3+fpvh*dconjg(funvh(-km,-kn))*spcom*sheff
    endif
22       continue
23      continue
      sum4hh=sum4hh+shh2*tempr1(n)*tempr3(m)
      sum4vv=sum4vv+svv2*tempr1(n)*tempr3(m)
      sum4hv=sum4hv+shv2*tempr1(n)*tempr3(m)
      sum4vh=sum4vh+svh2*tempr1(n)*tempr3(m)

      sum5hh=sum5hh+shh3*tempr2(n)*tempr2(m)
      sum5vv=sum5vv+svv3*tempr2(n)*tempr2(m)
      sum5hv=sum5hv+shv3*tempr2(n)*tempr2(m)
      sum5vh=sum5vh+svh3*tempr2(n)*tempr2(m)
!      if(cdabs(sum4vh-old).le.torlant) go to 30
!       old=sum4vh
24     continue
25    continue
30    continue
!---------------------------c
! sum up each term involved c
!---------------------------c
       sumhh=0.0d0
       sumvv=0.0d0
       sumhv=0.0d0
       sumvh=0.0d0
       rcom=0.125*k*k*dexp(-sig2*(ksz*ksz+kz*kz))
	   fphh=funhh(-kx,-ky)
       fpvv=funvv(-kx,-ky)
       fphv=funhv(-kx,-ky)
       fpvh=funvh(-kx,-ky)
       fshh=funhhs(-ksx,-ksy)
       fsvv=funvvs(-ksx,-ksy)
       fshv=funhvs(-ksx,-ksy)
       fsvh=funvhs(-ksx,-ksy) 
       sumhh=cdabs(fphh)**2*sum1+      &
			2.0*dreal(dconjg(fshh)*fphh)*sum2+     &
			cdabs(fshh)**2*sum3+(sum4hh+sum5hh)/pi
       coterm(1)=sumhh
       sumvv=cdabs(fpvv)**2*sum1+       &
            2.0*dreal(dconjg(fsvv)*fpvv)*sum2+     &
            cdabs(fsvv)**2*sum3+(sum4vv+sum5vv)/pi
       coterm(2)=sumvv
       sumhv=cdabs(fphv)**2*sum1+     &
             2.0*dreal(dconjg(fshv)*fphv)*sum2+     &
			 cdabs(fshv)**2*sum3+(sum4hv+sum5hv)/pi
       coterm(3)=sumhv
       sumvh=cdabs(fpvh)**2*sum1+      &
            2.0*dreal(dconjg(fsvh)*fpvh)*sum2+     &
	        cdabs(fsvh)**2*sum3+(sum4vh+sum5vh)/pi
       coterm(4)=sumvh 
35    continue
      do 40 ip=1,4
       sigma0(ip)=(katerm(ip)+crterm(ip)+coterm(ip))
40    continue
!------------------( ratio )-------------------------------------c
      if(ite.eq.3) then
         tranh=coterm(1)/sigma0(1)
         tranv=coterm(2)/sigma0(2)
      endif
!----------------------------------------------------------------c
      return
      end


!-------------------------------------------------------------------c
      function funvv(u,v)
      implicit real*8 (a-h,k,o-z)
      complex(8) k2,er,ur,funvv,rh,rv,ra
      complex(8) sq,sqs,Tv,Tvm
      common /cont1/ ksx,ksz,kx,kz,ksy,ky,si,sis,co,ccs, &
                    cp,sp,cps,sps,cpq,spq
      common /wlng/k,k2
      common /eu/er,ur
      common /rhhvv/rh,rv,ra

      s=si
      ss=sis
      cs=co
      css=ccs
      sf=spq
      csf=cpq
      sq=(ur*er-s*s)**0.5
      sqs=(ur*er-ss*ss)**0.5
      c1=(csf-s*ss)/(sq*css)
      c1s=(csf-s*ss)/(sqs*cs)
      c2=s*(ss-s*csf)/css
      c2s=ss*(s-ss*csf)/cs

      Tv=1+rv
      Tvm=1-rv

      funvv=-(cs*Tvm-sq*Tv/er)*(Tv*csf+Tvm*er*c1) &
          +(Tvm*Tvm-cs*Tv*Tvm/sq)*c2
      return
      end

      function funvvs(u,v)
      implicit real*8 (a-h,k,o-z)
      complex(8) k2,er,ur,funvvs,rh,rv,ra
      complex(8) sq,sqs,Tv,Tvm   
      common /cont1/ ksx,ksz,kx,kz,ksy,ky,si,sis,co,ccs, &
                    cp,sp,cps,sps,cpq,spq
      common /wlng/k,k2     
      common /eu/er,ur
      common /rhhvv/rh,rv,ra
       
      s=si
      ss=sis
      cs=co
      css=ccs
      sf=spq
      csf=cpq
      sq=(ur*er-s*s)**0.5
      sqs=(ur*er-ss*ss)**0.5
      c1=(csf-s*ss)/(sq*css) 
      c1s=(csf-s*ss)/(sqs*cs)
      c2=s*(ss-s*csf)/css
      c2s=ss*(s-ss*csf)/cs
         
      Tv=1+rv
      Tvm=1-rv
     
      funvvs=-(css*Tvm-sqs*Tv/er)*(Tv*csf+Tvm*er*c1s) &
           +(Tv*Tv-css*Tv*Tvm/sqs)*c2s
      return
      end

!------------------------------------------------------------------c
      function funhh(u,v)
      implicit real*8 (a-h,k,o-z)
      complex(8) k2,er,ur,funhh,rh,rv,ra
      complex(8) sq,sqs,Th,Thm
      common /cont1/ ksx,ksz,kx,kz,ksy,ky,si,sis,co,ccs, &
                    cp,sp,cps,sps,cpq,spq
      common /wlng/k,k2
      common /eu/er,ur
      common /rhhvv/rh,rv,ra

      s=si
      ss=sis
      cs=co
      css=ccs
      sf=spq
      csf=cpq
      sq=(ur*er-s*s)**0.5
      sqs=(ur*er-ss*ss)**0.5
      c1=(csf-s*ss)/(sq*css)
      c1s=(csf-s*ss)/(sqs*cs)
      c2=s*(ss-s*csf)/css
      c2s=ss*(s-ss*csf)/cs
         
      Th=1+rh
      Thm=1-rh

      funhh=(cs*Thm-sq*Th/ur)*(Th*csf+Thm*ur*c1) &
           -(Thm*Thm-cs*Th*Thm/sq)*c2
      return
      end


      function funhhs(u,v)
      implicit real*8 (a-h,k,o-z)
      complex(8) k2,er,ur,funhhs,rh,rv,ra
      complex(8) sq,sqs,Th,Thm
      common /cont1/ ksx,ksz,kx,kz,ksy,ky,si,sis,co,ccs, &
                    cp,sp,cps,sps,cpq,spq
      common /wlng/k,k2  
      common /eu/er,ur
      common /rhhvv/rh,rv,ra
      
      s=si
      ss=sis
      cs=co
      css=ccs
      sf=spq  
      csf=cpq
      sq=(ur*er-s*s)**0.5
      sqs=(ur*er-ss*ss)**0.5
      c1=(csf-s*ss)/(sq*css)
      c1s=(csf-s*ss)/(sqs*cs)
      c2=s*(ss-s*csf)/css
      c2s=ss*(s-ss*csf)/cs

      Th=1+rh
      Thm=1-rh
      
      funhhs=(css*Thm-sqs*Th/ur)*(Th*csf+Thm*ur*c1s)     &
	        -(Th*Th-css*Th*Thm/sqs)*c2s   
      return
      end


      function funvh(u,v)
      implicit real*8 (a-h,k,o-z)
      complex(8) k2,er,ur,funvh,rh,rv,ra
      complex(8) sq,sqs,Tp,Tm
      common /cont1/ ksx,ksz,kx,kz,ksy,ky,si,sis,co,ccs,     &
	                 cp,sp,cps,sps,cpq,spq
      common /wlng/k,k2
      common /eu/er,ur
      common /rhhvv/rh,rv,ra

      s=si
      ss=sis
      cs=co
      css=ccs 
      sf=spq  
      csf=cpq
      sq=(ur*er-s*s)**0.5
      sqs=(ur*er-ss*ss)**0.5
      c1=(csf-s*ss)/(sq*css)
      c1s=(csf-s*ss)/(sqs*cs)
      c2=s*(ss-s*csf)/css
      c2s=ss*(s-ss*csf)/cs
      
      Tp=1+ra
      Tm=1-ra

      funvh=(cs*Tp-sq*Tm/ur)*(Tm/css+Tp*ur/sq)*sf     &
	        + (Tp*Tp-cs*Tp*Tm/sq)*s*s*sf
      return
      end

      function funvhs(u,v)
      implicit real*8 (a-h,k,o-z)
      complex(8) k2,er,ur,funvhs,rh,rv,ra
      complex(8) sq,sqs,Tp,Tm
      common /cont1/ ksx,ksz,kx,kz,ksy,ky,si,sis,co,ccs,     &
	                 cp,sp,cps,sps,cpq,spq
      common /wlng/k,k2  
      common /eu/er,ur
      common /rhhvv/rh,rv,ra
      
      s=si
      ss=sis
      cs=co
      css=ccs
      sf=spq  
      csf=cpq
      sq=(ur*er-s*s)**0.5
      sqs=(ur*er-ss*ss)**0.5
      c1=(csf-s*ss)/(sq*css)
      c1s=(csf-s*ss)/(sqs*cs)
      c2=s*(ss-s*csf)/css
      c2s=ss*(s-ss*csf)/cs
      
      Tp=1+ra
      Tm=1-ra
      
      funvhs=-(css*Tm-sqs*Tp/er)*(Tp/cs+Tm*er/sqs)*sf     &
	        - (Tp*Tp-css*Tp*Tm/sqs)*ss*ss*sf  
      return
      end


      function funhv(u,v)
      implicit real*8 (a-h,k,o-z)
      complex(8) k2,er,ur,funhv,rh,rv,ra
      complex(8) sq,sqs,Tp,Tm
      common /cont1/ ksx,ksz,kx,kz,ksy,ky,si,sis,co,ccs,     &
	                 cp,sp,cps,sps,cpq,spq
      common /wlng/k,k2
      common /eu/er,ur
      common /rhhvv/rh,rv,ra

      s=si
      ss=sis
      cs=co
      css=ccs
      sf=spq
      csf=cpq
      sq=(ur*er-s*s)**0.5   
      sqs=(ur*er-ss*ss)**0.5
      c1=(csf-s*ss)/(sq*css)
      c1s=(csf-s*ss)/(sqs*cs)
      c2=s*(ss-s*csf)/css
      c2s=ss*(s-ss*csf)/cs
         
      Tp=1+ra
      Tm=1-ra

      funhv=(cs*Tm-sq*Tp/er)*(Tp/css+Tm*er/sq)*sf     &
	        + (Tm*Tm-cs*Tp*Tm/sq)*s*s*sf
      return
      end

      function funhvs(u,v)
      implicit real*8 (a-h,k,o-z)
      complex(8) k2,er,ur,funhvs,rh,rv,ra
      complex(8) sq,sqs,Tp,Tm
      common /cont1/ ksx,ksz,kx,kz,ksy,ky,si,sis,co,ccs,     &
	                 cp,sp,cps,sps,cpq,spq
      common /wlng/k,k2  
      common /eu/er,ur
      common /rhhvv/rh,rv,ra
      
      s=si
      ss=sis
      cs=co
      css=ccs
      sf=spq 
      csf=cpq
      sq=(ur*er-s*s)**0.5
      sqs=(ur*er-ss*ss)**0.5
      c1=(csf-s*ss)/(sq*css)
      c1s=(csf-s*ss)/(sqs*cs)
      c2=s*(ss-s*csf)/css
      c2s=ss*(s-ss*csf)/cs
      
      Tp=1+ra
      Tm=1-ra
      
      funhvs=-(css*Tp-sqs*Tm/ur)*(Tm/cs+Tp*ur/sqs)*sf     &
	        - (Tm*Tm-css*Tp*Tm/sqs)*ss*ss*sf  
      return
      end

!-------------------------------------------------------------c


!-----------------------------------c
!  surface roughness spectrum       c
!-----------------------------------c
      function spectrum(u,v,fn)
      implicit real*8 (a-h,k,o-z)
      dimension zt(512),wt(512)
      external funtexp