SUBROUTINE Search_center(data_in,nx,ny,rx,ry, * xc_o,yc_o,xc_n,yc_n,max_val,IFIND,TDO) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C This subroutine uses the downhill scheme to search for C the maximun centre of data_in. C this subroutine maybe used for searching vortex centre,... C C C INPUT: C data_in Input data C nx,ny Dimensions of all matrices and vectors. C rx,ry Coordinate vectors in rad. C xc_o,yc_o First guess centre. C TDO Half centre search domain C C OUTPUT C xc_n,yc_n New vorticity centre in rad. C max_val Maximum value at center C IFIND Flag for centre detection (= 0 if no distinct C vorticity extremum was found). C C C NOTE: unit used: meter C if you want to use with other unit, please modify HDOWN, HDOWNMIN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMPLICIT COMPLEX (Z) C C REAL data_in(nx,ny),Coef(nx,ny,4,4) REAL rx(nx),ry(ny) REAL WW(3) real max_val C COMPLEX AA,VV,UU(7),Z(3) C C C...PARAMETERS AND DETERMINATION OF THE SEARCH DOMAIN C DX=Rx(2)-Rx(1) PP=0.1 ! Must be greater than -1 HDOWN=2.*Dx ! First step size downhill method (deg) HDHMIN=0.1*Dx ! Minimum step size downhill method (deg) NFIN = 100 ! Max. no. of downhill iterations WMIN=-1.E20 ! Boundary value downhill search IFIND=1 TDOR=TDO XU=xc_o-TDOR XO=xc_o+TDOR YU=xc_o-TDOR YO=yc_o+TDOR ! test ! write(2,*)'data_in, in search_center' ! do j=1,ny ! write(2,'(150f10.2)')(data_in(i,j),i=1,nx) ! enddo ! write(2,*)'rx',rx ! write(2,*)'xc_o',xc_o ! write(2,*)'ry',ry ! write(2,*)'yc_o',yc_o C C...COMPUTE THE COEFFICIENT MATRIX FOR BIRATIONAL INTERPOLATION C IERR=0 CALL RAT2D(data_in,rx,ry,nx,ny,nx,ny,PP,Coef,IERR) IF (IERR.EQ.1) THEN PRINT*,' Search: NO CENTRE SEARCH POSSIBLE' IFIND=0 RETURN ENDIF C C...INITIAL SUBSTITUTIONS FOR DOWNHILL METHOD C ZS=CMPLX(xc_o,yc_o) UU(1)=(1.,0.) UU(2)=(.8660254,.5) UU(3)=(0.,1.) UU(4)=(.9659258,.2588190) UU(5)=(.7071068,.7071068) UU(6)=(.2588190,.9659258) UU(7)=(-.2588190,.9659258) HIN=HDOWN NEND=0 C C...CALCULATION OF MIDPOINT VALUE C CALL RATVA(Coef,rx,ry,nx,ny,nx,ny,PP, & REAL(ZS),AIMAG(ZS),WO) C C...SEARCH FOR MAXIMUM C 201 KK=1 II=0 202 VV=(-1.,0.) 203 AA=(-.5,.8660254) 204 Z(1)=ZS+HIN*VV*AA DXZ=REAL(Z(1))-xc_o DYZ=AIMAG(Z(1))-yc_o RZ=SQRT(DXZ*DXZ+DYZ*DYZ) IF (RZ.GT.TDOR) THEN WW(1)=WMIN ELSE CALL RATVA(Coef,rx,ry,nx,ny,nx,ny,PP, & REAL(Z(1)),AIMAG(Z(1)),WW(1)) ENDIF Z(2)=ZS+HIN*VV DXZ=REAL(Z(2))-xc_o DYZ=AIMAG(Z(2))-yc_o RZ=SQRT(DXZ*DXZ+DYZ*DYZ) IF (RZ.GT.TDOR) THEN WW(2)=WMIN ELSE CALL RATVA(Coef,rx,ry,nx,ny,nx,ny,PP, & REAL(Z(2)),AIMAG(Z(2)),WW(2)) ENDIF Z(3)=ZS+HIN*CONJG(AA)*VV DXZ=REAL(Z(3))-xc_o DYZ=AIMAG(Z(3))-yc_o RZ=SQRT(DXZ*DXZ+DYZ*DYZ) IF (RZ.GT.TDOR) THEN WW(3)=WMIN ELSE CALL RATVA(Coef,rx,ry,nx,ny,nx,ny,PP, & REAL(Z(3)),AIMAG(Z(3)),WW(3)) ENDIF NEND=NEND+1 IF (NEND.GT.NFIN) THEN !PRINT*,' CENTR: Iteration number exceeded' IFIND=0 IF (WO-WW(1)) 224,223,223 223 IF (WO-WW(2)) 227,225,225 225 IF (WO-WW(3)) 222,219,219 224 IF (WW(1)-WW(2)) 227,227,226 226 IF (WW(1)-WW(3)) 222,220,220 227 IF (WW(2)-WW(3)) 222,221,221 219 GOTO 118 220 ZS=Z(1) GOTO 118 221 ZS=Z(2) GOTO 118 222 ZS=Z(3) GOTO 118 ENDIF IF (WW(1)-WW(3)) 206,205,205 205 IF (WW(1)-WW(2)) 208,208,207 206 IF (WW(2)-WW(3)) 209,208,208 207 NNR=1 GOTO 210 208 NNR=2 GOTO 210 209 NNR=3 210 IF (WO-WW(NNR)) 212,212,211 211 GOTO (213,214,215),KK 212 KK=1 II=0 AA=(.7071068,.7071068) VV=(Z(NNR)-ZS)/HIN WO=WW(NNR) ! write(*,*)"WO=",WO ! write(*,*)"HIN=",HIN ! write(*,*)"NEND=",NEND ZS=Z(NNR) GOTO 204 213 KK=2 IF (HIN.LT.HDHMIN) GOTO 118 HIN=HIN*.25 GOTO 203 214 KK=3 HIN=HIN*4. GOTO 202 215 II=II+1 IF (II-7) 216,216,217 216 VV=UU(II) GOTO 203 217 IF (HIN.LT.HDHMIN) GOTO 118 HIN=HIN*.25 II=0 GOTO 202 118 xc_n=REAL(ZS) yc_n=AIMAG(ZS) max_val=WO C c bh2 ! write(*,'(" So buoc lap :",i5)')NEND ! write(*,'(" Do chinh xac:",f5.2,"km")') HIN*110 ! write(*,'(" Lon:",f6.2,", Lat:",f6.2)') xc_n,yc_n ! write(*,'(" Gia tri:",f7.2)') -WO RETURN END