-
Notifications
You must be signed in to change notification settings - Fork 0
/
antiKsitogravfd.f90
82 lines (82 loc) · 3.32 KB
/
antiKsitogravfd.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
subroutine AntiKsitogravfd(calcpntfl,boundryfl,ksigrdfl,dr)
implicit none
character*800::calcpntfl,boundryfl,ksigrdfl
character*800::line,str,astr
integer sn,len,astat(8),i,j,nlon,nlat,kk
real*8::dr,hd(6),hd1(6),hd2(6),rec(800),pi,RAD,mr,rr
real*8::GRS(6),BLH(3),BLH0(3),rln(3),NFD(5),gr
real*8::PoissonBLH,rst,db,vv,sk0,P1,P2,gra,rga,vms,vmw
real*8,allocatable::ksi(:,:),sfh(:,:)
integer::status=0
!---------------------------------------------------------------------
pi=datan(1.d0)*4.d0;RAD=pi/180.d0;mr=36.d2/RAD
GRS(1)= 3.986004415d14; GRS(2)=6378136.3d0; GRS(3)=1.082636277388d-3
GRS(4) = 7.292115d-5; GRS(5) = 1.d0/298.25641153d0
open(unit=8,file=ksigrdfl,status="old",iostat=status)
if(status/=0)goto 902
read(8,'(a)') line
call PickRecord(line,len,rec,sn)
hd(1:6)=rec(1:6)
nlat=nint((hd(4)-hd(3))/hd(6))
nlon=nint((hd(2)-hd(1))/hd(5))
hd(5)=(hd(2)-hd(1))/real(nlon)
hd(6)=(hd(4)-hd(3))/real(nlat)
allocate(ksi(nlat,nlon), stat=astat(1))
allocate(sfh(nlat,nlon), stat=astat(2))
if (sum(astat(1:2)) /= 0) then
close(8);goto 902
endif
do i=1,nlat
read(8,*,end=903)(ksi(i,j),j=1,nlon)
enddo
903 close(8)
open(unit=10,file=boundryfl,status="old",iostat=status)
if(status/=0)goto 904
read(10,'(a)') line
call PickRecord(line,len,rec,sn)
hd1(1:6)=rec(1:6)
if(sum(hd1-hd)>1.d-5)then !格网规格不同The grid specifications are different
close(10);goto 904
endif
do i=1,nlat
read(10,*,end=905)(sfh(i,j),j=1,nlon)
enddo
905 close(10)
db=hd(5)/16.d0;vv=hd(5)*1.2d4
open(unit=8,file=calcpntfl,status="old",iostat=status)
if(status/=0)goto 904
open(unit=10,file='reslt.txt',status="replace")
read(8,'(a)') line !读取头文件read the file header
write(10,101)trim(line)
kk=0
do while(.not.eof(8))
read(8,'(a)') line
call PickRecord(line,len,rec,sn)
if(sn<4)goto 906
kk=kk+1;BLH(2)=rec(2);BLH(1)=rec(3);BLH(3)=rec(4)!the ellipsoidal height
call BLH_RLAT(GRS,BLH,rln);call GNormalfd(BLH,NFD,GRS);gr=NFD(2);rr=rln(1)
sk0=PoissonBLH(BLH,ksi,sfh,nlat,nlon,hd,dr,GRS)
BLH0=BLH;BLH0(3)=BLH(3)+vv/2.d0
P1=PoissonBLH(BLH0,ksi,sfh,nlat,nlon,hd,dr,GRS)
BLH0=BLH;BLH0(3)=BLH(3)-vv/2.d0
P2=PoissonBLH(BLH0,ksi,sfh,nlat,nlon,hd,dr,GRS)
rga=(P2-P1)/vv*gr*1.d5;gra=rga-gr*sk0/rr/2.d0*1.d5
BLH0=BLH;BLH0(1)=BLH(1)+db/2.d0
P1=PoissonBLH(BLH0,ksi,sfh,nlat,nlon,hd,dr,GRS)
BLH0=BLH;BLH0(1)=BLH(1)-db/2.d0
P2=PoissonBLH(BLH0,ksi,sfh,nlat,nlon,hd,dr,GRS)
vms=(P2-P1)/(rr*db*RAD)*mr
BLH0=BLH;BLH0(2)=BLH(2)+db/dcos(rln(2)*RAD)/2.d0
P1=PoissonBLH(BLH0,ksi,sfh,nlat,nlon,hd,dr,GRS)
BLH0=BLH;BLH0(2)=BLH(2)-db/dcos(rln(2)*RAD)/2.d0
P2=PoissonBLH(BLH0,ksi,sfh,nlat,nlon,hd,dr,GRS)
vmw=(P2-P1)/(rr*db*RAD)*mr
write(10,101)trim(line),gra,rga,vms,vmw
if(kk/200*200==kk)write(*, '(a,i9)'), ' Calculated point number: ',kk
enddo
906 close(8);close(10)
904 deallocate(ksi,sfh)
902 continue
write (*,*)' Complete the computation! The results are saved in the file reslt.txt.'
101 format(a,40F12.4)
end