sub_indat_p.f
4.96 KB
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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
c======================================================================|
subroutine sub_indat(fnin,fntmp,angref6,nang,dtr,nnset
1 ,ifnpo,idprop)
c======================================================================|
implicit none
integer,intent(in) :: nang,ifnpo,idprop
real*8,intent(in) :: angref6(nang),dtr
character*100 ,intent(in):: fntmp(nang)
integer,intent(out) :: nnset
integer :: i,j,ifn0,ifn1,ifn2,iang
real*8 :: parain(12),js1
real*8 :: angref1,dtin,dtset
character*100 :: fnin,fnout
integer :: nnin
real*8,allocatable :: datain(:,:),jsin(:),dataset(:,:),jsset(:)
real*8,allocatable :: dataset1(:)
real*8,allocatable :: dang(:),jstmp(:)
real*8,allocatable :: dang2(:),dang3(:),jstmp2(:),jstmp3(:)
real*8 :: dangtmp0,dtmp1,paramx(8)
ifn0=ifnpo
ifn1=41
ifn2=42
!---count line number
call count_lines(fnin,nnin)
allocate(datain(13,nnin),jsin(nnin),dang(nnin),jstmp(nnin))
allocate(dang2(nnin),dang3(nnin),jstmp2(nnin),jstmp3(nnin))
!---read input data
open(ifn1,file=fnin,status='old')
do i=1,nnin
read(ifn1,*) js1,(parain(j),j=1,12)
!!!!!!! elena : Dimensions mismatch (13/12)
datain(1:12,i)=parain(:)
jsin(i)=js1
enddo
dtin=(jsin(2)-jsin(1))*float(idprop) !--REV0318
!---set output time
dtset=dtr
nnset=int(float(nnin)*dtin/dtset)
allocate(dataset(13,nnset),jsset(nnset),dataset1(nnset))
do i=1,nnset
jsset(i)=jsin(1)+dtset*float(i)*float(idprop) !--REV0318
enddo
!---shift to 6 reference longitude & output
do iang=1,nang !---6 ref.angles
angref1=angref6(iang)
dang=datain(10,:)-angref1
!---angle interpolation (note +/-180 deg.)
dang2(1)=dang(1)
do i=2,nnin
dangtmp0=dang(i)-dang(i-1)
if(abs(dangtmp0).ge.300.) dangtmp0=0.
dang2(i)=dang2(i-1)+dangtmp0
enddo
dang3=dang2
if(sum(dang2)/float(nnin).gt.360.) dang2=dang2-360.
if(sum(dang2)/float(nnin).gt.360.) dang2=dang2-360.
if(sum(dang2)/float(nnin).lt.-360.) dang2=dang2+360.
if(sum(dang2)/float(nnin).lt.-360.) dang2=dang2+360.
if(sum(dang2)/float(nnin).lt.0.) dang3=dang2+360.
if(sum(dang2)/float(nnin).ge.0.) dang3=dang2-360.
do i=1,nnin
if(dang(i).lt.-180.) dang(i)=dang(i)+360.
if(dang(i).gt.180.) dang(i)=dang(i)-360.
enddo
jstmp=jsin-dang/360.*26.*86400.
jstmp2=-dang2/360.*26.*86400.
jstmp3=-dang3/360.*26.*86400.
call convim(jstmp,jsset,jstmp2,dataset1
1 ,nnin,nnset,jstmp2(1),jstmp2(nnin))
dataset(2,:)=dataset1
call convim(jstmp,jsset,jstmp3,dataset1
1 ,nnin,nnset,jstmp3(1),jstmp3(nnin))
dataset(3,:)=dataset1
do i=1,nnset
j=2
if(abs(dataset(2,i)).gt.abs(dataset(3,i))) j=3
dataset(13,i)=dataset(j,i)
enddo
c---angle jump (1/2)
dang=datain(10,:)
dang2(1)=dang(1)
do i=2,nnin
dangtmp0=dang(i)-dang(i-1)
if(abs(dangtmp0).ge.300.) dangtmp0=0.
dang2(i)=dang2(i-1)+dangtmp0
enddo
datain(10,:)=dang2
dang=datain(12,:)
dang2(1)=dang(1)
do i=2,nnin
dangtmp0=dang(i)-dang(i-1)
if(abs(dangtmp0).ge.300.) dangtmp0=0.
dang2(i)=dang2(i-1)+dangtmp0
enddo
datain(12,:)=dang2
c---data conversion
do i=1,12
call convim(jstmp,jsset,datain(i,:),dataset1
1 ,nnin,nnset,datain(i,1),datain(i,nnin))
dataset(i,:)=dataset1
c---angle jump (2/2)
if((i.eq.10.or.i.eq.12))then
do j=1,nnset
if(dataset(i,j).lt.0.) dataset(i,j)=dataset(i,j)+360.
if(dataset(i,j).lt.0.) dataset(i,j)=dataset(i,j)+360.
if(dataset(i,j).gt.360.) dataset(i,j)=dataset(i,j)-360.
if(dataset(i,j).gt.360.) dataset(i,j)=dataset(i,j)-360.
enddo
endif
enddo
c---interpolation using non-extreme value
do i=2,nnset
dangtmp0=dataset(10,i)-dataset(10,i-1)
dataset1(i)=dangtmp0
enddo
dataset1(1)=dataset1(2)
dtmp1=sum(abs(dataset1))/float(nnset)
do j=1,8
paramx(j)=sum(abs(dataset(j,:)))/float(nnset)+1.
enddo
do i=1,nnset
dangtmp0=dataset1(i)
if(abs(dangtmp0).le.dtmp1*0.01)then
do j=1,8
if(abs(dataset(j,i)).gt.abs(paramx(j))*5.)
1 dataset(j,i)=paramx(j)
enddo
endif
enddo
c---data output
fnout=fntmp(iang)
open(ifn2,file=fnout,status='unknown',form='formatted')
do i=1,nnset
write(ifn2,'(f12.0,13(1pe15.5))') jsset(i),dataset(1:13,i)
enddo
close(ifn2,status='keep')
enddo !---6 ref.angles
c---closing
close(ifn1,status='keep')
deallocate(jsin,jsset,datain,dataset,dang,jstmp)
return
end
c==================================================