cxdef 1440 linear -179.875 0.25 cydef 720 linear -89.875 0.25 ! include "../../include/all.inc" 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) 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) character *7 tail 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 (10,file="dsource_is.asc",form="unformatted",status="old") read (10) s k=0 do i=1,ill do j=1,jll if(s(i,j).gt.5.)s(i,j)=0. enddo enddo hots=maxval(s) c-------------------------------------------------------------- c open (10,file="dsource_is.asc",status="old") c read (10,*) s close(10) print *,minval(s),maxval(s),im,jm mask=1 UND = -99 imah=0 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).eq.hots) 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 endif else if(s(i1,j1).le.0.1) 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 endif else if(s(i1,j1).ge.0.19.and.s(i1,j1).le.0.21) then call hbox ( ar7,IM,JM,olak,olok,i2,j2) if(i2.ne.0.and.j2.ne.0) then ima3(i2,j2)=ima3(i2,j2)+1 endif else if(s(i1,j1).ge.0.39.and.s(i1,j1).le.0.41) then call hbox ( ar7,IM,JM,olak,olok,i2,j2) if(i2.ne.0.and.j2.ne.0) then ima4(i2,j2)=ima4(i2,j2)+1 endif else if(s(i1,j1).ge.0.64.and.s(i1,j1).le.0.66) then call hbox ( ar7,IM,JM,olak,olok,i2,j2) if(i2.ne.0.and.j2.ne.0) then ima5(i2,j2)=ima5(i2,j2)+1 endif endif enddo enddo sumh= sum(SUM(imah, DIM=1),dim=1) sum3= sum(SUM(ima3, DIM=1),dim=1) sum4= sum(SUM(ima4, DIM=1),dim=1) sum5= sum(SUM(ima5, DIM=1),dim=1) sum0= sum(SUM(ima0, DIM=1),dim=1) sum100= sum(SUM(ima100, DIM=1),dim=1) print *,'sum0,sum3,sum4,sum5,sumh,sum100' print *,sum0,sum3,sum4,sum5,sumh,sum100 print *,'maxval ima,ima3,ima4,ima5,>2 imah,ima100' print *,maxval(ima0),maxval(ima3),maxval(ima4),maxval(ima5) print *,maxval(imah),maxval(ima100) prag=maxval(ima100)/2. a3=0.2 a4=0.4 a5=0.65 c ah=1.8*(a3*sum3+a4*sum4+a5*sum5)/sumh ah=(a3*sum3+a4*sum4+a5*sum5)/sumh do i=1,im do j=1,jm if(ima100(i,j).ge.prag)then src(i,j) =0. else write(88,*)"--",i,j,ima100(i,j),imah(i,j) write(88,*)ima3(i,j),ima4(i,j),ima5(i,j) src(i,j)=a3*ima3(i,j)+a4*ima4(i,j)+a5*ima5(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 enddo enddo c mxsrc=maxval(src) mxsrc=ah*maxval(ima0) print *,"ah=",ah,"maxval=",mxsrc do j=1,jm do i=1,im s2(i,j)=src(i,j)/mxsrc ! if(s2(i,j).le.0.006)s2(i,j)=100. if(s2(i,j).le.0.006)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 print*, "minval maxval s2 ", minval(s2(:,:)), maxval(s2(:,:)) rmax=maxval(s2(:,:)) print*, "rmax rmax rmax ", rmax s2(:,:)=s2(:,:)/rmax where(s2(:,:).lt.0.1)s2(:,:)=s2(:,:)*1.6 s2(:,:)=s2(:,:)*3. print*, "FINAL maska minval maxval s2 ", + minval(s2(:,:)), maxval(s2(:,:)) !feb2019predisland open(20,file="../../output/PREF.DAT",form="unformatted") open(20,file="PREF.DAT",form="unformatted") write (20) s2 close (20) 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) 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