bri_tem_nonsnow.f90

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

subroutine BRI_TEM_NONSNOW(Ang,Fre,Co_deg,Co_len,Sv,Cv,Fv,&
                           
			   Dz,Tss,Wliq,Wice,Tve,SinScaAlb,b,Lai,Tb)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!the BT procedure is used to simulate the brightness temperature 
!received by the passive microwave sensors. 
!---------------------------------------------------------------
!output parameters:
!  Tb: simulation brightness temperatures received by sensors
!      including H polarization Tb and V polarization Tb at two
!      different frencencies
!---------------------------------------------------------------
!input parameters:
!  Ang: observation angle [degrees]
!  Fre: observation freqencies [GHz]
!  Co_deg: coarseness degree [cm]
!  Co_len: correlation length [cm]
!  Sv: sand volume percentage
!  Cv: clay volume percentage
!  Fv: solid materials volume percentage
!  Dz: snow and soil layer depth in CLM [m]
!  Tss:temperature in each layer for snow and soil in CLM [K]
!  Wliq: water content in each layer for snow and soil in CLM [kg/m2]
!  Wice: ice   content in each layer for snow and soil in CLM [kg/m2]
!  Tve:  vegetation effective temperature [K]
!  SinScaAlb: single scattering albedo
!  b: vegetation parameter used to determine optical path for some freqency
!     for more dails about this parameter,refer to articles by Jackson
!  Lai: a parameter for vegetation used to parameterize Wv [kg/m2]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none

real,parameter::ConCoe=0.035  ! conversion coefficent from Lai to Wv, unit[kg/m2] 

real,dimension(1:2)::Tb  ! output paramter brightness temperatures for 
                         ! H and V polarizations,Tb(1) vertical,Tb(2) horizontal[k]

real::Ang,Fre,Co_deg,Co_len,Sv,Cv,Fv,Tve,SinScaAlb,b,Lai

real,dimension(-4:10)::Dz,Tss,Wliq,Wice ! inputting into EFF_VAL_COM for computing Tb 

real::&     ! mediate variables outputting from EFF_VAL_COM  
      Tse,& ! Tse --soil effective temperature
      Vms,& ! Vms --soil effective volumetric water content           
      Vis   ! Vis --soil effective volumetric ice   content		  
				  

real::Trans	! vegetation canopy transmittance for micowave at some frequency

real,dimension(1:2)::Emi ! the H and V emittances for bare soil 
                         ! at some frequency
						 ! Emi(1) for vertical,Emi(2) for horizontal
						 ! derive from IEM

real::Tspv,&  ! bare soil vertically polarized brightness temperature
      Tsph    ! bare soil horizontally polarized brightness temperature
       
!print*, Ang 
!print*, fre 
!print*, co_deg 
!print*, co_len 
!print*, sv 
!print*, cv 
!print*, Sv
!print*, fv 
!print*, dz 
!print*, Tss
!print*, wliq 
!print*, wice
!print*, tve
!print*,sinscaalb
!print*, b
!print*, lai
!print*, tb 

!print*,'Tse',Tse,'Vms',vms,'Vis',Vis
call EFF_VAL_COM( Dz,Tss,Wliq,Wice,Ang,Fre,Sv,Cv,Fv,Tse,Vms,Vis)
!print*,'Tse',Tse,'Vms',vms,'Vis',Vis
!print*,'ang', ang
call IEM(Ang,Fre,Co_deg,Co_len,Sv,Cv,Fv,Tse,Vms,Vis,Emi)
!print*,'Emi',Emi
!print*,'ang', ang

Tspv=Emi(1)*Tse
Tsph=Emi(2)*Tse

Trans=exp(-b*(Lai*ConCoe)/cos(Ang))

Tb(1)=(1+(1-Emi(1))*Trans)*(1-Trans)*(1-SinScaAlb)*Tve+ &
      Trans*Tspv! +5.0

Tb(2)=(1+(1-Emi(2))*Trans)*(1-Trans)*(1-SinScaAlb)*Tve+ &
      Trans*Tsph! +5.0


continue


! the detailed description of above formula could be found in the article
! by Jackon in 1991 and 1993

end subroutine BRI_TEM_NONSNOW
!========================================================================


subroutine EFF_VAL_COM( Dz,Tss,Wliq,Wice,Ang,Fre,Sv,Cv,Fv,Tse,Vms,Vis)

!      ------------------------------------------------------------------
!     |                                                                  |
!     |  EFF_VAL_COM is used to compute effive temperature,volume water  |
!     |  content and volume ice content in conjunction with CLM. For the |
!     |  time being,only top 4 soil layers are employed in this          |
!     |  procedure. However, more soil layers can be added to it. For    |
!     |  detailed formula, please refer to atbd-amsr-land.pdf by Njouk   |
!     |   the article by Jianchen Shi  and CLM by Yongjiu Dai            |
!     |                                                                  |
!     |  first version time :5/20/2004   author: QinJun                  |
!     |  1th modified version time: 7/23/2004  author: QinJun            |
!     |                                                                  |
!     |                                                                  |
!     |                                                                  |
!     |                                                                  |
!      ------------------------------------------------------------------
implicit none

real,parameter:: Rhow= 1.00*10**3   ! water density [kg/m3]

real,parameter:: Rhoi= 0.91*10**3   ! ice density [kg/m3]

real,parameter:: Pi  = 3.1415926  

integer,parameter:: Ln=5               ! the number of layers used to compute effective values

real,dimension(-4:10):: Dz, &  !Dz: snow and soil layer depth in CLM [m]
                        Tss,&  !temperature in each layer for snow and soil in CLM [K]
                        Wliq,& !water content in each layer for snow and soil in CLM [kg/m2]
                        Wice   !ice content in each layer for snow and soil in CLM [kg/m2]

real:: Ang,&    ! sensor zenith angle
       Fre      ! sensor frequency

real:: Sv,& ! Sv -- sand volume percentage
       Cv,& ! Cv -- clay volume percentage
	   Fv   ! Fv -- percentage of solid soil excluding ice

real:: Tse,& ! soil effective temperature
       Vms,& ! soil effective volume moisture content
       Vis   ! soil effective volume ice      content

real,dimension(1:Ln)::Alpha ! parameters used for computing effective values

real:: Ext_Coe   ! extinction coefficient used for comuting effective values

real:: Epsilonr,&   ! real part in soil dielectric constant 
       Epsiloni     ! imagory part in soil dielectric constant

real,dimension(1:Ln):: Wvfrac,&     ! temporary varaible denoting volume water content for each soil layer in CLM
                       Ivfrac       ! temporary varaible denoting volume ice   content for each soil layer in CLM

real,dimension(1:300):: STss,   &   ! sublayer temperature used to computing Tse
                        SWvfrac,&
					    SIvfrac,&
					    SAlpha ,&
					    SDz    

                    
real:: iv

real:: wn ! wavenumber

integer:: i,j,k,l,m,n     ! iteration cycle number

Tse=0.0E+0;Vms = 0.0E+0; Vis = 0.0E+0 

wn=(2*Pi)*Fre/0.3    ! transforming Frequency to wavenumber

do i=1,Ln

!   if (i.EQ.1) then

      Wvfrac(i)=(Wliq(i)/Rhow)/Dz(i)

      Ivfrac(i)=(Wice(i)/Rhoi)/Dz(i)

!   else

!      Wvfrac(i)=((Wliq(i)-Wliq(i-1))/Rhow)/Dz(i) ! converting CLM layer water content to Epsion volume water content

!      Ivfrac(i)=((Wice(i)-Wice(i-1))/Rhoi)/Dz(i) ! converting CLM layer ice   content to Epsion volume ice   content

!   end if 

   call epsion(Fre,Sv,Cv,Fv,Tss(i),Wvfrac(i),Ivfrac(i),Epsilonr,Epsiloni,iv)

   Alpha(i)=2 * wn * &
            abs(aimag(csqrt(dcmplx(Epsilonr,Epsiloni)-(sin(Ang))**2)))

end do


do m=1,300

    if (m.LE.50) then

       STss(m)=Tss(1)

	   SWvfrac(m)=Wvfrac(1)

	   SIvfrac(m)=Ivfrac(1)

	   SAlpha (m)=Alpha (1)

	   SDz    (m)=Dz(1)/50.0

	 else if (m.GT.50.AND.m.LE.100) then

           STss(m)=Tss(2)

	   SWvfrac(m)=Wvfrac(2)

	   SIvfrac(m)=Ivfrac(2)

	   SAlpha (m)=Alpha (2)

	   SDz    (m)=Dz(2)/50.0	 
	  
	 else if (m.GT.100.AND.m.LE.180) then

           STss(m)=Tss(3)

	   SWvfrac(m)=Wvfrac(3)

	   SIvfrac(m)=Ivfrac(3)

	   SAlpha (m)=Alpha (3)

	   SDz    (m)=Dz(3)/80.0	
	  
	 else if  (m.GT.180.AND.m.LE.240) then	
	    
	   STss(m)=Tss(4)

	   SWvfrac(m)=Wvfrac(4)

	   SIvfrac(m)=Ivfrac(4)

	   SAlpha (m)=Alpha (4)

	   SDz    (m)=Dz(4)/60.0

	 else if  (m.GT.240.AND.m.LE.300) then

	   STss(m)=Tss(5)

	   SWvfrac(m)=Wvfrac(5)

	   SIvfrac(m)=Ivfrac(5)

	   SAlpha (m)=Alpha (5)

	   SDz    (m)=Dz(5)/60.0

	 !else if  (m.GT.1820.AND.m.LE.4620) then

	   !STss(m)=Tss(6)

	   !SWvfrac(m)=Wvfrac(6)

	   !SIvfrac(m)=Ivfrac(6)

	   !SAlpha (m)=Alpha (6)

	   !SDz    (m)=Dz(6)/2800.0

	! else if  (m.GT.4620.AND.m.LE.6930) then

	  ! STss(m)=Tss(7)

	  ! SWvfrac(m)=Wvfrac(7)

	   !SIvfrac(m)=Ivfrac(7)

	  ! SAlpha (m)=Alpha (7)

	  ! SDz    (m)=Dz(7)/2310.0	   
	   
	 end if 	
	      
end do



do j=1,300
      
    Ext_Coe=0.0
      
    do k=1,j
      
       Ext_Coe=Ext_Coe+SAlpha(k)*SDz(k)
	     
    end do

      Tse=Tse+STss (j)*SAlpha(j)*exp(-Ext_Coe)*SDz(j)  
      
      Vms=Vms+SWvfrac(j)*SAlpha(j)*exp(-Ext_Coe)*SDz(j)

      Vis=Vis+SIvfrac(j)*SAlpha(j)*exp(-Ext_Coe)*SDz(j)

end do 

return

end subroutine EFF_VAL_COM




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine IEM(ang,fre,co_deg,co_len,sv,cv,fv,tse,vms,vis,emi) 
! be used to couple with CLM
! ff is an array used to store silumation frenencies in GHz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!c
!                                                           c
!  program computes bistatics scattering coefficients       c
!  from "3d" finitely conducting surface                    c
!  ( only for HH and VV polarization )                      c
!  ------------------------------------------------------   c
!  input parameters:                                        c
!     er: surface relative dielectric constant              c
!     cl: surface correlation length in cm                  c
!    sig: surface rms height in cm                          c
!     fr: operating frequency in GHz                        c
!  -----------------------------------------------------    c
!    itype: select what type of surface correlation         c
!           =1  Gaussian correlation function               c
!           =2  exponential correlation function            c
!           =3  transformed exponential correlation         c
!           =4 generalized power law spectrum               c
!              p = 1.5 exponential                          c
!              p = infinity  Gaussina                       c
!    icross:                                                c
!           =0 no cross terms are computed                  c
!           =1    corss terms are computed                  c
!    icort:                                                 c
!           =0 no corrected terms are computed              c
!           =1    corrected terms are computed              c
!    muti:                                                  c
!           =0 no mutiple scattering are computed           c
!           =1    mutiple scattering are computed           c 
!    -------------------------------------------------      c
!    npts: number of quadature points for integration       c
!    npol: number of polarizations to be computed           c
!    -------------------------------------------------      c
!    Arguments for subroutine "sigma" (on input):           c
!    ang: incident angle in deg                             c
!    angs: scattering angle in deg                          c
!    phi: incident azimuth angle in deg                     c
!    phis: scattering azimuth angle in deg                  c
!    z1: zeros in k direction for Gaussian quaduture        c
!    w1: weights in k direction for Gaussian quadature      c
!    z2: zeros in phi direction for Gaussian quadature      c
!    w2: weights in phi direction for Gaussian quadature    c
! rmslp: surface rms slope                                  c
!   ------------------------------------------------------  c
!    Arguments for subroutine "sigma" (on output):          c
!    sigma: array for IEM scattering coefficients           c
!          sigma(1) ---> HH polarization                    c
!          sigma(2) ---> VV polarization                    c
!          sigma(3) ---> HV polarization                    c
!          sigma(4) ---> VH polarization                    c
!    -----------------------------------------------------  c
!    ite=1:Fresnel reflection coeff. approxmated by R(ang)  c
!    ite=2:Fresnel reflection coeff. approxmated by R(0)    c
!    ite=3:calculate the proposed transition function       c
!    ite=4:Fresnel reflection coeff. approxmated by R(T)    c
!    -----------------------------------------------------  c
!    tranh:the ratio of the complementary term to the total c
!          scattering coefficients in HH polarization       c
!    tranv:the ratio of the complementary term to the total c
!          scattering coefficients in VV polarization       c
!    tranh0,tranv0: when k approaches zero                  c
!    -----------------------------------------------------  c
!    trfachh,trfacvv: proposed transition function          c
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!c 
!
!             Comment by ZHaoKaiguang                       !
!   when k*sig grows large, some variables will overflow,   !
!such as  tempk,tempr1 and etc. So in the program,only these!
!variables' LOG value is stored into the corresopding varia-!
!les whose name is the original ones puls "_LOG". To ensure !
!it to run properly, other parts are aslo modified.         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!      implicit real*8 (a-h,k,o-z)
       implicit real*8 (a-h,k,o-z)
      complex(8) er,ur,k2,ffhh,ffvv,stem,rh2,rv2, &
               gghh,ggvv,funhh,funvv,funvvs,funhhs
      complex(8) rhi,rvi
      dimension sigma0(4),sigma1(4)
      dimension w1(512),z1(512),w2(512),z2(512)
      dimension wt(512),zt(512)
      dimension shad_cor(512,512)
      logical  back,shadowing
      common /qua/zt,wt,npts
      common /cont2/ rxx,sig,sig2
      common /wlng/k,k2
      common /p/pi,pi2,spi,spi2
      common /al/cl,effslop,npol
      common /type/itype
      common /eu/er,ur
      common /index/icross,icort,muti
      integer npol,npts,ite
      real*8 uth, cu, fu
      common /order/uth,cu,fu
      real:: vms
  !    real*8 ang
!define variable names for epsion
      !real,dimension(2)::fre	! simulation frenencies
      real::fre
      real,dimension(1:2)::emi	! emittance of two bands first vertical second horizontal
      real::sv,cv,fv,tse,co_deg,co_len	 !for detailed specifications, see epsion 
      integer::ii        !iteration cycle number for frenencies
!print*,ang,fre,co_deg,co_len,sv,cv,fv,tse,vms,vis

!print*,'icross',icross,'icort',icort,'muti',muti
!
! -------------------------------------
! ** assign initial values to common **
! -------------------------------------

  npts = 0

  rxx = 0.0E+0;sig = 0.0E+0;sig2 = 0.0E+0

  k = 0.0E+0; k2 = (0.0E+0,0.0E+0)

  pi = 0.0E+0; pi2 = 0.0E+0; spi = 0.0E+0; spi2 = 0.0E+0

  cl = 0.0E+0; effslop = 0.0E+0; npol = 0

  itype = 0

  er = (0.0E+0,0.0E+0); ur = (0.0E+0,0.0E+0)

  icross = 0; icort = 0; muti = 0

   




!	  st1=st-273.16  ! convert from unit K to unit degree C for 
	                ! the requirement of Epsion model use to compute 
					! dielectic constant for bare soil

!!!!!!!!      open(15,file='iem.in',status='old',access='sequential')
      shadowing=.false.
      back=.false.
      pi=dacos(-1.0d0)
      pi2=2.0*pi
      ur=1.0
      npol=4			 !!polarization numbers hh,vv,hv,vh
      npts=512
      spi=sqrt(pi)
      spi2=sqrt(pi2)

!-------------------------------------------------------c
!         Input Parameters                              c
!-------------------------------------------------------c
      icross=1	!!cross terms are computed
      icort=1	!!correct terms are computed?
      muti=0	!!single scattering
      itype=1	!!Gaussian correlation function
      ite=2		!!Fresnel reflection coeff. approxmated by R(ang)