      subroutine HGradientTwo(dt,hgt,nlat,nlon,hd,rt1,rt2,CRS)
      !格网格值2阶水平梯度估计,x轴指向北,y轴指向东
      implicit none
	integer nlon,nlat,i,j,ki,kj,si,sj
	real*8::hd(6),dt(nlat,nlon),hgt(nlat,nlon),rt1(nlat,nlon),rt2(nlat,nlon)
	real*8::RAD,pi,rr,rlat,lon,rlat0,lon0,dlat,dlon,dx,dy
	real*8::sin2f,cos2f,sinf,cosa,sina
	real*8::BLH(3),rln(3),BPB(5,5),BB(5),BPL(5),xx(5),obs,dl,tmp
	real*8::CRS(5)
!---------------------------------------------------------------------
   	pi=datan(1.d0)*4.d0;RAD=pi/180.d0
      rt1=0.d0;rt2=0.d0
 	do i=1,nlat
        BLH(1)=hd(3)+(real(i)-0.5d0)*hd(6)
        do j=1,nlon
          BLH(2)=hd(1)+(real(j)-0.5d0)*hd(5)
          BLH(3)=hgt(i,j)
          call BLH_RLAT(CRS,BLH,rln)
	    rlat0=rln(2)*RAD;lon0=rln(3)*RAD
	    rr=rln(1);BPB=0.d0;BPL=0.d0;xx=0.d0
	    do ki=i-1,i+1
	      BLH(1)=hd(3)+(real(ki)-0.5d0)*hd(6)
	      do kj=j-1,j+1
	        BLH(2)=hd(1)+(real(kj)-0.5d0)*hd(5)
	        if(ki<1.or.ki>nlat.or.kj<1.or.kj>nlon)goto 9200
	        if(ki==i.and.kj==j)goto 9200
		    BLH(3)=hgt(ki,kj)
              call BLH_RLAT(CRS,BLH,rln)
	        rlat=rln(2)*RAD;lon=rln(3)*RAD
	        dlon=(lon-lon0)/2.d0;dlat=(rlat-rlat0)/2.d0
	        sin2f=dsqrt(dsin(dlat)**2+dsin(dlon)**2*dcos(rlat0)*dcos(rlat))
	        dl=2.d0*rr*sin2f;cos2f=dsqrt(1.d0-sin2f**2);sinf=2.d0*sin2f*cos2f
	        cosa=(dcos(rlat0)*dsin(rlat)-dsin(rlat0)
     >		    *dcos(rlat)*dcos(dlon*2.d0))/sinf
	        sina=dcos(rlat)*dsin(dlon*2.d0)/sinf
	        BB(1)=dl*cosa;BB(2)=dl*sina;obs=dt(ki,kj)-dt(i,j)
	        BB(3)=BB(1)*BB(2);BB(4)=BB(1)**2/2.d0;BB(5)=BB(2)**2/2.d0
	        do si=1,5
	          BPL(si)=BPL(si)+BB(si)*obs/dl
	          do sj=1,5
	            BPB(si,sj)=BPB(si,sj)+BB(si)*BB(sj)/dl
	          enddo
	        enddo
9200	        continue
	      enddo
	    enddo
	    tmp=0.d0
	    do si=1,5
	       do sj=1,si-1
	         BPB(sj,si)=BPB(si,sj)
	       enddo
             tmp=tmp+BPB(si,si)**2
	    enddo
	    tmp=dsqrt(tmp/5.d0)
	    do si=1,5
	      BPB(si,si)=BPB(si,si)+tmp*1.d-5
	    enddo
	    call EqHestens(BPB,xx,5,BPL)
          rt1(i,j)=xx(4);rt2(i,j)=xx(5)
        enddo
      enddo
      end
