cxdef 1440 linear -179.875 0.25
cydef 720 linear -89.875 0.25
include "all.inc"
PARAMETER(IM=-WBD/DLMD+1.5,JM=-2.*SBD/DPHD+1.5)
c parameter(ILL=1440, JLL=720)
c parameter(bw= -179.875, BS=-89.875,DLL=0.25)
C parameter(ILL=7200, JLL=3600)
C parameter(bw= -179.975, BS=-89.975,DLL=0.05)
parameter(ILL=14400, JLL=7200)
parameter(bw= -179.9875, BS=-89.9875,DLL=0.025)
parameter(imh=9,ip=8)
dimension id7(7),ar7(7), s(ILL,JLL),s2(im,jm)
dimension imah(im,jm),ima0(im,jm),src(im,jm),ima100(im,jm)
dimension ima3(im,jm),ima4(im,jm),ima5(im,jm),mask(im,jm)
dimension dist(im,jm,ip),imah8(im,jm,imh),gpss(imh,ip)
dimension ihot(im,jm),rhot(im,jm),sums5(im,jm),imas5(im,jm)
character *7 tail
character *1 dum1
i0w=5000
i0e=7700
j0s=5050
j0n=7200
c hots=1.65
DTR = .01745329
AR7(1)=2. !0 - latlon; 1 - imjm; 2 - im,jm
AR7(2)=TLM0D
AR7(3)=WBD
AR7(4)=DLMD
AR7(5)=TPH0D
AR7(6)=SBD
AR7(7)=DPHD
ctph0=cos(tph0d*dtr)
stph0=sin(tph0d*dtr)
c-----------island
open (11,file="datad.csv",form="formatted",status="old")
read(11,*)
do i=1,imh
read(11,*) dum1,dum1,(gpss(i,j),j=1,ip)
print *, "****",dum1,dum1,(gpss(i,j),j=1,ip)
enddo
close(11)
open (10,file="dsource_s5_is.asc",form="unformatted",status="old")
read (10) s
print *,"minval=",minval(s)," maxval=",maxval(s)
k=0
do i=1,ill
do j=1,jll
if(s(i,j).gt.115.)s(i,j)=0.
enddo
enddo
hots=6 !maxval(s)
c--------------------------------------------------------------
c open (10,file="dsource_is.asc",status="old")
c read (10,*) s
close(10)
c print *,minval(s),maxval(s),im,jm
mask=1
UND = -99
imah=0
imah8=0
ihot(:,:)=9
do i1=i0w,i0e
olok=bw+(i1-1)*dll
do j1=j0s,j0n
olak=bs+(j1-1)*dll
call hbox ( ar7,IM,JM,olak,olok,i2,j2)
if(i2.ne.0.and.j2.ne.0) then
ima0(i2,j2)=ima0(i2,j2)+1
endif
if(s(i1,j1).ge.6..and.s(i1,j1).le.13.) then
call hbox ( ar7,IM,JM,olak,olok,i2,j2)
if(i2.ne.0.and.j2.ne.0) then
imah(i2,j2)=imah(i2,j2)+1
is8=int(s(i1,j1))-5 !hotspot je 6,7,8,9,10,11,12,13
imah8(i2,j2,is8)= imah8(i2,j2,is8)+1
ihot(i2,j2)=is8
endif
! else if(s(i1,j1).le.0.0001) then
else if(s(i1,j1).ge.17.) then
call hbox ( ar7,IM,JM,olak,olok,i2,j2)
if(i2.ne.0.and.j2.ne.0) then
ima100(i2,j2)=ima100(i2,j2)+1
write(77,*),i2,j2,ima100(i2,j2)
endif
else if(s(i1,j1).gt.0.0001.and.s(i1,j1).le.1.) then
call hbox ( ar7,IM,JM,olak,olok,i2,j2)
if(i2.ne.0.and.j2.ne.0) then
imas5(i2,j2)=imas5(i2,j2)+1
sums5(i2,j2)=sums5(i2,j2)+s(i1,j1)
endif
endif
enddo
enddo
! write(33,*)maxval(ima0),maxval(imas5),maxval(ima100)
! &,maxval(imah)
! do i=30,54
! do j=85,110
! write(33,*)ima0(i,j),imas5(i,j),ima100(i,j),imah(i,j)
! write(33,*)sums5(i,j),imah(i,j)
! end do
! end do
! stop
do i=1,im
do j=1,jm
ih=ihot(i,j)
if(ih.ne.9)print *,"========",i,j,ih,(gpss(ih,l),l=1,8)
do k=1,ip
dist(i,j,k)=gpss(ih,k)
if(imah8(i,j,k).ne.0.and.imah(i,j).ne.imah8(i,j,k)) then
print *,"*PAZNJA* tacki ",i,j," pripada vise hots "
& , imah8(i,j,k)," iz hotspota ",k
endif
enddo
enddo
enddo
open (20,file="particule_is.bin",form="unformatted")
write(20)dist
close(20)
sumh= sum(SUM(imah, DIM=1),dim=1)
sum3= sum(SUM(imas5, DIM=1),dim=1)
sumas5= sum(SUM(sums5, DIM=1),dim=1)
sum0= sum(SUM(ima0, DIM=1),dim=1)
sum100= sum(SUM(ima100, DIM=1),dim=1)
print *,'sum0,sum3,sumh,sum100'
print *,sum0,sum3,sumh,sum100
print *,'maxval ima,ima3,>2 imah,ima100'
print *,maxval(ima0),maxval(imas5)
print *,maxval(imah),maxval(ima100)
prag=maxval(ima100)/10.
c ah=1.8*(a3*sum3+a4*sum4+a5*sum5)/sumh
ah=sumas5/(sumh*4.)
print *,"ah=",ah
maxis5=maxval(imas5)
mxsrc=ah*maxval(imah)
do i=1,im
do j=1,jm
if(ima100(i,j).ge.prag)then
write(33,*)ima0(i,j),imas5(i,j),ima100(i,j)
src(i,j) =0.
else
c write(88,*)
src(i,j)=sums5(i,j)+ah*imah(i,j)
endif
c if(imah(i,j).ne.0..or.ima3(i,j).ne.0..or.ima4(i,j).ne.0..or.
c & ima5(i,j).ne.0.)then
c write(55,*)i,j,imah(i,j),ima3(i,j),ima4(i,j),ima5(i,j)
c write(55,*)"==",imah(i,j)+ima3(i,j)+ima4(i,j)+ima5(i,j),ima0(i,j)
c endif
write(88,*)ima0(i,j),ima100(i,j)
&,imah(i,j),imas5(i,j),src(i,j)
enddo
enddo
c mxsrc=maxval(src)
c mxsrc=ah*maxval(ima0)
rmv=maxval(src)!/mxsrc
print *,"ah=",ah,"maxval=",mxsrc ,"rmv=",rmv
do j=1,jm
do i=1,im
s2(i,j)=src(i,j)!/mxsrc
s2(i,j)=s2(i,j)/rmv
if(s2(i,j).le.0.005)s2(i,j)=0.
c write(87,*)ima3(i,j),ima4(i,j),ima5(i,j),imah(i,j)
write(87,*)s2(i,j)
tph = sbd+(j-1)*dphd
tlm = wbd+mod(j+1,2)*dlmd+(i-1)*2*dlmd
tlm = tlm * dtr
tph = tph * dtr
call RTLLMx(TLM,TPH,TLM0D,DTR,CTPH0,STPH0,ro,ra)
c
p = (ro-BW)/DLL
q = (ra-BS)/DLL
i0 = p + 1
j0 = q + 1
p = p - int(p)
q = q - int(q)
z1 = s (i0 ,j0)
z2 = s (i0+1,j0)
z3 = s (i0+1,j0+1)
z4 = s (i0 ,j0+1)
c vall2(ii)=blx(p,q, z1, z2, z3, z4)
ccc if(imah(i,j).ne.0) then
c print *,"=========",imah(i,j),ima0(i,j)
ccc s2(i,j)=1.*imah(i,j)/ima0(i,j)
ccc print *,i,j,"==1==",s2(i,j)
ccc else
ccc s2(i,j)=blu (UND,p,q,z1, z2, z3, z4)/(1.*ima0(i,j))
ccc print *,i,j,"==0==",s2(i,j)
ccc endif
c print *,s2(i,j),(hots*ima0(i,j)),imah(i,j)
enddo
enddo
c
open (20,file="PREF.DAT",form="unformatted")
write (20) s2
close (20)
D=-99.
ID7(1)=11
ID7(2)=01
ID7(3)=15
ID7(4)=00
ID7(5)=0
ID7(6)=00 ! minf
ID7(7)=00 ! secf
print *,ID7
print *,"_________________",im,jm
print *,AR7
print *,"____________",maxval(s2)
CALL WGRADS(ID7,36,0,0,1,0,AR7,IM,JM,1,0.,S2 ,DUM)
rhot=ihot/1.
CALL WGRADS(ID7,35,0,0,1,0,AR7,IM,JM,1,0.,rhot ,DUM)
do l=1,ip
lev = l
print *,l,lev,"part=",maxval(dist(:,:,l))
CALL WGRADS(id7,37,0,0,100,lev,ar7,im,jm,1,1.*lev,dist(:,:,l),D)
enddo
stop
c
open (10,file="tails.txt",status="old")
dxy = 1./120.
1 read(10,"(a)",end=99) tail
read (tail(2:4),*) bwt
read (tail(6:7),*) bnt
if(tail(1:1).eq."W") bwt = -bwt
if(tail(5:5).eq."S") bnt = -bnt
jj = 0
do rlat=bnt,bnt-50,-dxy
jj = jj + 1
ii = 0
do rlon=bwt,bwt+40,+dxy
ii = ii + 1
enddo
enddo
print *,tail,"xx",jj,ii
goto 1
99 close(10)
end
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
SUBROUTINE RTLLMx(TLM,TPH,TLM0D,DTR,CTPH0,STPH0,ALMD,APHD)
C
STPH=SIN(TPH)
CTPH=COS(TPH)
CTLM=COS(TLM)
STLM=SIN(TLM)
C
APH=ASIN(STPH0*CTPH*CTLM+CTPH0*STPH)
CPH=COS(APH)
C
ALMD=TLM0D+ASIN(STLM*CTPH/CPH)/DTR
APHD=APH/DTR
C
RETURN
END