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