bri_tem_nonsnow.f90

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

!!!!!!!!	  501   read(15,*,end=1998) fr,ang,cl,sig,err,eri,vms



	!do 1998 ang=    55.0,		 56,	  50	  ! ang: observation zenith anlge
	!do 1998 ff=     10.65,		11.44,	  50      ! frequency
	!do 1998 vms=     0.3,		 0.5,    10.02    ! volume soil moisture
	!do 1998 sig=     0.25,		 3.0,	 0.25	  ! sigma: rms Height
	!do 1998 cl=      5,			20.0,	  55	  ! Correlation length
	sig=co_deg ! coarseness degree
	cl =co_len ! correlation length
	do ii=1,    1,   1

	call epsion(fre,sv,cv,fv, tse,vms,vis,err,eri,vi)
	!fre:frenency,sv:sand,cv:clay,fv:bulk density,st:soil temperature&
	!vms:volume soil moisture,vis:volume soil ice,err:real part,erri:imaginary part
	!vi:none

!-------------------------------------------------------c 
      k=pi2*fre/30.0         !! wave number: wave length unit is cm
      kl=k*cl		        !! horizontal roughness presentation
      ksig=k*sig	        !! Vertical roughness presentation
      er=dcmplx(err,eri)	!! complex type of dielectric constant
      cl2=cl*cl
      k2=k*cdsqrt(er*ur)	
      sig2=sig*sig
      if(itype.eq.1) slope=1.414*sig/cl		!!
      if(itype.eq.2) slope=2.0*sig/cl
      if(itype.eq.3) slope=1.732*sig/cl
      rmslp=slope
      slope=rmslp
      effslop=0.0+slope
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! "effslop" is artifical rms slope used to take
! second-order shadowing into account
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! if mutiple scattering is desired, then
! generate gaussian quad. for integration
!
      if(muti.eq.0) go to 16	! go to 16 for single scattering
!------------------------------------------c
!    set up integration limits             c
!------------------------------------------c
       a1=0.0
       b1=(1.0-1.0e-20)*k
       a2=0.0
       b2=pi
!-----------------------------------------c
!  generates zeros and absissa            c
!-----------------------------------------c
      call quagen(zt,wt,npts)		 ! subroutine is at line 1112
      as1=(b1-a1)/2.0
      bs1=(b1+a1)/2.0
      as2=(b2-a2)/2.0
      bs2=(b2+a2)/2.0
      do 10 i=1,npts
       w1(i)=wt(i)*as1
       z1(i)=zt(i)*as1+bs1
       w2(i)=wt(i)*as2
       z2(i)=zt(i)*as2+bs2
 10   continue
!-------------------------------------------c
! compute and store 2nd-order shadowing     c
! function,then pass to subroutine          c
!-------------------------------------------c
       do 15 i=1,npts
       do 15 j=1,npts
		u=z1(i)*dcos(z2(j))
		v=z1(i)*dsin(z2(j))
! 2nd-order incident angle determined by "u" and "v"
		ql=dasind(dsqrt(u**2+v**2)/k)
		call shadow(back,ql,ql,shef)
		shad_cor(i,j)=shef*shef
15     continue   
16    continue				 ! the above is for Multiscattering

!---------------------------------------------------------------c
!              main program                                     c
!---------------------------------------------------------------c


      co=dcosd(ang)

      si=dsind(ang)

      phi=0.0
      angs = ang
      phis = 179.99

      ehh=0.0
      evv=0.0
      ehv=0.0
      evh=0.0
      ehht=0.0
      evvt=0.0
      evht=0.0
      ehvt=0.0

      do 20 is = 1,89,2		! step for incidence angle
         angs = is
      do 20 js = 1,179,2	! step for azimuth angle
         phis = js
      if(phis.eq.180.0) phis=179.99		 ! ??
      if(phis.eq.90.0) phis=89.99		 ! ??
      if(angs.eq.0.0) angs=0.01			 ! ??

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

      call sigma(trfachh,trfacvv,ang,angs,phi,phis,z1,w1,z2,w2,	&
      npts,shad_cor,sigma0,tranh,tranv,rmslp,ite,rhi,rvi)



!      call sigmaT(trfachh,trfacvv,ang,angs,phi,phis,z1,w1,z2,w2,&
!      npts,shad_cor,sigma1,rmslp,ite,rhi,rvi)

!----------------------------------------------------------------c
       if(.not.shadowing) then
         shfct=1.0
       else
         call shadow(back,abs(ang),abs(angs),shfct)
       endif
       sh2=shfct**2

      if(sigma0(1).gt.0.0) hh=(sigma0(1)*sh2)
      if(sigma0(2).gt.0.0) vv=(sigma0(2)*sh2)
      if(sigma0(3).gt.0.0) hv=(sigma0(3)*sh2)
      if(sigma0(4).gt.0.0) vh=(sigma0(4)*sh2)

      if(sigma1(1).gt.0.0) hht=(sigma1(1)*sh2)
      if(sigma1(2).gt.0.0) vvt=(sigma1(2)*sh2)
      if(sigma1(3).gt.0.0) hvt=(sigma1(3)*sh2)
      if(sigma1(4).gt.0.0) vht=(sigma1(4)*sh2)

!----------------------------------------------------------------c
!     calculate transition functions in backscattering           c
!----------------------------------------------------------------c
       if(ite.eq.3) then      ! pass because ite=1       
       stem=cdsqrt(er)
       rv2=(er-stem)/(er+stem)			 ! Fresnel coe. of horizontal pol.
       rh2=(1.0-stem)/(1.0+stem)		 ! Frensel coe. of Vertical pol.
	   ffhh=-2.0*rh2/co
       ffvv=2.0*rv2/co
       vvv=0.0
	   gghh=(funhh(-k*si,vvv)+funhhs(k*si,vvv))		! compute the Fpp
       ggvv=(funvv(-k*si,vvv)+funvvs(k*si,vvv))
      tranh0=											  &
           0.25*cdabs(gghh)**2.0/((cdabs(ffhh)**2.0)*4.0  &
           +(dconjg(ffhh)*gghh+ffhh*dconjg(gghh))+		  &
           0.25*cdabs(gghh)**2.0)					 ! ?
      tranv0= &
           0.25*cdabs(ggvv)**2.0/((cdabs(ffvv)**2.0)*4.0  &
           +(dconjg(ffvv)*ggvv+ffvv*dconjg(ggvv))+        &
           0.25*cdabs(ggvv)**2.0)					 ! ?


      trfachh=(1.0-tranh/tranh0)
      trfacvv=(1.0-tranv/tranv0)
      if(trfachh.gt.1) trfachh=1.0					 ! ?
      if(trfachh.lt.0) trfachh=0.0					 ! ?
      if(trfacvv.gt.1) trfacvv=1.0					 ! ?
      if(trfacvv.lt.0) trfacvv=0.0					 ! ?

      endif
!----------------------------------------------------------------c
      ehh=ehh+abs(hh*dsind(angs)/(4.0*pi*dcosd(ang)))		  ! for ite not equal to 3
      evv=evv+abs(vv*dsind(angs)/(4.0*pi*dcosd(ang)))
      evh=evh+abs(vh*dsind(angs)/(4.0*pi*dcosd(ang)))
      ehv=ehv+abs(hv*dsind(angs)/(4.0*pi*dcosd(ang)))

      ehht=ehht+abs(hht*dsind(angs)/(4.0*pi*dcosd(ang)))
      evht=evht+abs(vht*dsind(angs)/(4.0*pi*dcosd(ang)))
      evvt=evvt+abs(vvt*dsind(angs)/(4.0*pi*dcosd(ang)))
      ehvt=ehvt+abs(hvt*dsind(angs)/(4.0*pi*dcosd(ang)))

20    continue

       if(.not.shadowing) then
       shfct=1.0
       else
        call shadow(back,abs(ang),abs(ang),shfct)
      endif

       sh2=shfct**2
       rh0=cdabs(rhi)**2
       rv0=cdabs(rvi)**2
       cohh=rh0*dexp(-(2*k*dcosd(ang)*sig)**2)*shfct**2
       covv=rv0*dexp(-(2*k*dcosd(ang)*sig)**2)*shfct**2

       covvt=dexp(-sig2*real(  &
       ((k2*k2-k*k+(k*dcosd(ang))**2)**0.5-k*dcosd(ang))**2.0))

       cohht=(1-rh0)*covvt*shfct**2
       covvt=(1-rv0)*covvt*shfct**2

      evv=evv*(pi/90.0)**2.0
      ehv=ehv*(pi/90.0)**2.0
      ehh=ehh*(pi/90.0)**2.0
      evh=evh*(pi/90.0)**2.0

      evvt=evvt*(pi/90.0)**2.0*CDABS(cdsqrt(er))
      ehvt=ehvt*(pi/90.0)**2.0*CDABS(cdsqrt(er))
      ehht=ehht*(pi/90.0)**2.0*CDABS(cdsqrt(er))
      evht=evht*(pi/90.0)**2.0*CDABS(cdsqrt(er))

      evt=1-(2*evv+2*ehv+covv)
      eht=1-(2*ehh+2*evh+cohh)

      evtt=(2*evvt+2*ehvt+covvt)
      ehtt=(2*ehht+2*evht+cohht)

! added by zhaokg      
      cch=rh0*dexp(-(2*k*dcosd(ang)*sig)**2)
      ccv=rv0*dexp(-(2*k*dcosd(ang)*sig)**2)
! added by zhaokg

	  !open(unit=3,file='IEMOUT.TXT',status='unknown',&
	  !    access='sequential',form='formatted')
      !write(3,99) fr, ang, vms, sig,cl,evt,eht	   !evt:v plorization eht: h plorization
!	  write(3,99) 1-evt,1-eht,1-evtt/((1-evt)+evtt),1-ehtt/((1-eht)+ehtt), sig, cl, rv0, rh0,vms,covv,cohh,cv,ch
      emi(1)=evt  
	  emi(2)=eht

99    format(14f10.6)

117   format(2(f8.6,1x),2(f8.6,1x),2(f5.2,1x),f6.3,1x,f4.1,1x,&
            2(f7.4,1x))
!!!!!      go to 501

!1998  continue
end do


end subroutine IEM
	  !contains

!*******************************************************************c
! subroutine sigma calculates the scattering coefficients           c
!*******************************************************************c
      subroutine sigma(trfachh,trfacvv,sita,sitas,phi,phis,z1,w1,z2, &
                      w2,nw,shad_cor,sigma0,tranh,tranv,rmslp,ite,	 &
                      rhi,rvi)
      implicit real*8 (a-h,k,o-z)
      integer nw,ite
      complex(8) k2,er,ur,rv,rh,ra,ssem,siem,stem,fvv,fhh,fvh,fhv
      complex(8) rh0,rv0,fkh,fkv,rhi,rvi
      complex(8) shh1,shh2,shh3,svv1,svv2,svv3
      complex(8) shv1,shv2,shv3,svh1,svh2,svh3
      complex(8) sum3hh,sum3vv,sum3hv,sum3vh
      complex(8) sum4hh,sum4vv,sum4hv,sum4vh
      complex(8) sum5hh,sum5vv,sum5hv,sum5vh
      complex(8) funhh,funvv,funhv,funvh
      complex(8) funhhs,funvvs,funhvs,funvhs
      complex(8) fphh,fpvv,fphv,fpvh
      complex(8) fshh,fsvv,fshv,fsvh
      complex(8) old
      dimension katerm(4),crterm(4),coterm(4),sigma0(4)
      dimension tempk(2),tempc1(2),tempc2(2)
	  dimension tempk_LOG(30000),tempc1_LOG(30000),tempc2_LOG(30000)
      dimension tempr1(1),tempr2(1),tempr3(1)
	  dimension tempr1_LOG(30000),tempr2_LOG(30000),tempr3_LOG(30000)
	  dimension z1(512),w1(512),z2(512),w2(512)
      dimension shad_cor(512,512)
      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 /type/itype
      common /rhhvv/rh,rv,ra
      common /cont2/ rxx,sig,sig2
      common /p/pi,pi2,spi,spi2
      common /al/cl,effslop,npol
      common /index/ icross,icort,muti
!~~~~~~~~~~~~~~~~~~~~~~~~~~c
! clear all arrays         c
!~~~~~~~~~~~~~~~~~~~~~~~~~~c
      kl=k*cl
!      do 5 i=1,15000			 ! the number should be kept uniform with dimension
!       tempk(i)=0.0
!       tempr1(i)=0.0
!       tempr2(i)=0.0
!       tempr3(i)=0.0
!       tempc1(i)=0.0
!       tempc2(i)=0.0
!5     continue
      torlant=1.0e-18
      co=dcosd(sita)
      ccs=dcosd(sitas)
      si=dsind(sita)
      si2=si*si
      sis=dsind(sitas)
!      if(sis.lt.0.0) return
      cp=dcosd(phi)
      cps=dcosd(phis)
      sp=dsind(phi)
      sps=dsind(phis)
      cpq=dcosd(phis-phi)
      spq=dsind(phis-phi)
!---------------------------c
      kz=k*co
      ksz=k*ccs
      kx=k*si*cp
      ky=k*si*sp
      ksx=k*sis*cps
      ksy=k*sis*sps
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
! reflection coefficients are based on the incid angles c
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
      stem=cdsqrt(er*ur-si*si)
      ssem=cdsqrt(er*ur-sis*sis)
      siem=cdsqrt(er*ur-si2)
      rvi=(er*co-stem)/(er*co+stem)
      rhi=(ur*co-stem)/(ur*co+stem)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
! rv0,rh0 are evaluated at normal incidence             c
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
      rv0=(er-cdsqrt(er*ur))/(er+cdsqrt(er*ur))
      rh0=(ur-cdsqrt(er*ur))/(ur+cdsqrt(er*ur))
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~c
      if(ite.eq.1) then
        rh=rhi
        rv=rvi
      endif
      if(ite.eq.2) then
        rh=rh0
        rv=rv0
      endif
      if(ite.eq.3) then
        rh=rh0
        rv=rv0
!       new and need to be deleted
! rh=rhi
! rv=rvi
      endif
      if(ite.eq.4) then
        rh=rhi+(rh0-rhi)*trfachh
        rv=rvi+(rv0-rvi)*trfacvv
      endif
      ra=(rv-rh)/2.0
!---------------------------------------------c
!     kirchhoff field coefficients            c
!---------------------------------------------c
      f=(si*sis-(1.0+co*ccs)*cpq)/(ccs+co)
      fvv=2.0*rv*f
      fhh=-2.0*rh*f
      fhv=2.0*ra*spq
      fvh=-2.0*ra*spq
      fkh=-2.0*rh0*f
      fkv=2.0*rv0*f
!-------------------------c
! compute kirchhoff term  c
!-------------------------c
      tempk(1)= (sig*(ksz+kz))**2

	  tempk_LOG(1)=dlog(tempk(1))

      n=2
	  do  2 n=2,1.5* tempk(1)+100
      
1     fn=float(n)
!      tempk(n)=tempk(1)*tempk(n-1)/fn
!      if(dabs(tempk(n)-tempk(n-1)).le.torlant) go to 2
!      n=n+1  
!      go to 1  
  
  	   tempk_LOG(n)=tempk_LOG(1)+tempk_LOG(n-1)-dlog(fn)
2     continue
!      mm=n
       mm=n-1
!      print*,mm,tempk(mm)
!      write(*,*) mm
      sumk=0.0
	  sumk_LOG=0
	  
	  ep_LOG=-sig2*(ksz+kz)**2+dlog(0.5*k*k)

      do 3 n=1,mm
!	  sumk=sumk+tempk(n)*spectrum(ksx-kx,ksy-ky,fn)

	  fn=float(n)
	  sumk_LOG=sumk_LOG+dexp(tempk_LOG(n)+ep_LOG)*spectrum(ksx-kx,ksy-ky,fn)

3     continue  

!     ep=0.5*k*k*exp(-sig2*(ksz+kz)**2)
!     sumkep=sumk*ep

      sumkep=sumk_LOG 
	  fn=1
      katerm(1)=sumkep*cdabs(fhh)**2
      katerm(2)=sumkep*cdabs(fvv)**2
      katerm(3)=sumkep*cdabs(fhv)**2
      katerm(4)=sumkep*cdabs(fvh)**2
      krh=sumkep*cdabs(fkh)**2
      krv=sumkep*cdabs(fkv)**2

!------------------------------------------c
!  end of kirchhoff term computation       c
!------------------------------------------c

      if(ite.eq.3) then
!       if(ite.gt.1) then
        rh=rh0
        rv=rv0
      else
        rh=rhi
        rv=rvi
      endif
!      ra=(rv-rh)/2.0

      if(icross.eq.0) go to 17
!------------------------------------------c
! compute cross term (with 3 sub-terms)    c
!------------------------------------------c
!	crcom=0.5*k*k*exp(-sig2*(ksz*ksz+kz*kz+ksz*kz))
    crcom_LOG=-sig2*(ksz*ksz+kz*kz+ksz*kz)+dlog(0.5*k*k)
	tempc1(1)=sig2*ksz*(ksz+kz)
	tempc2(1)=sig2*kz*(ksz+kz)
	tempc1_LOG(1)=dlog(tempc1(1))
	tempc2_LOG(1)=dlog(tempc2(1))
	do 10 n=2,mm
	fn=float(n)
! 	tempc1(n)=tempc1(1)*tempc1(n-1)/fn
    tmp_LOG=dlog(fn)
  	tempc1_LOG(n)=tempc1_LOG(1)+tempc1_LOG(n-1)-tmp_LOG
	tempc2_LOG(n)=tempc2_LOG(1)+tempc2_LOG(n-1)-tmp_LOG
10   continue   
!		do 11 n=2,mm
!		 fn=float(n)
!		 tempc2(n)=tempc2(1)*tempc2(n-1)/fn
!11		continue
	
	sum1=0.0
	sum2=0.0

	sum1_LOG=0.0
	sum2_LOG=0.0

	sum3vv=(0.0,0.0)
	sum3hh=(0.0,0.0)
	sum3hv=(0.0,0.0)
	sum3vh=(0.0,0.0)
	old=(0.0,0.0)
	do 12 n=1,mm
	fn=float(n)
	spe=spectrum(ksx-kx,ksy-ky,fn)
!	sum1=sum1+tempc1(n)*spe
!	sum2=sum2+tempc2(n)*spe
!	tmp_LOG=dlog(spe)
	sum1_LOG=sum1_LOG+dexp(tempc1_LOG(n)+crcom_LOG)*spe
	sum2_LOG=sum2_LOG+dexp(tempc2_LOG(n)+crcom_LOG)*spe
12    continue

	sum1=sum1_LOG
	sum2=sum2_LOG

if(muti.eq.0) go to 16
!----------------------------------------c
! calculate the mutiple scattering term  c
!----------------------------------------c
do 15 n=1,mm
   fn=float(n)
do 14 m=1,mm
    fm=float(m)
  shh1=(0.0,0.0)
  svv1=(0.0,0.0)
  shv1=(0.0,0.0)
  svh1=(0.0,0.0)
  do 13 i=1,nw
   do 13 j=1,nw
    u=z1(i)*dcos(z2(j))
    v=z1(i)*dsin(z2(j))
     ql=dasind(dsqrt(u**2+v**2)/k)
!        call shadow(back,ql,ql,sheff)