! ! Topical Cyclone Detection Program for regional climate models ! (C) Bui Hoang Hai, 2010 ! Trinh Tuan Long,2012 ! Kim Nguyen, 2012 use datedef Implicit none !Dimension defines Integer :: Nx, Ny Integer, parameter :: Np=4 Real :: plev(Np) Data plev /850.,700.,500.,300./ Real, allocatable, dimension(:) :: Rx, Ry Real :: Rx0, Dx character(len=20):: Ryfile,input_dir ! Variables on Pressure levels and sea level, ptp: potential temp Real, allocatable, dimension(:,:,:) :: up,vp,wp,tp,hp,Vorp,ptp Real, allocatable, dimension(:,:) :: Slp1, Slp2, Vtang ! Variables for criteria and threholds Real, dimension(Np) :: Tanop Real :: Tcenter,Hcenter,Pcenter,Tout,Pout,Hout,SLPano,Tano,TCoreano,h500ano,OCS, Relax ! These one will be passed from the namelist real :: VorCrit, & ! Vorticicy criteria, e.g. 1.e-4 [s-1] TCorecrit, & ! Core Temp.anomaly on , e.g. > 1K H500crit, & ! Geopoential Heigh Ano at 500 mb, e.g < -10 OCScrit, & ! OCS of 850hpa criteria, e.g. 5 [m/s] SLPcrit, & ! Sea Level anomaly, e.g. < -3 hPa CPradius, & ! Radius for detecting the same vortices (current and the pass timeslice), e.g. 3.5 degree Cradius ! Radius for detecting the same vortices (current time slice), e.g. 3.5 degree ! Pressure center Searching variables Real, parameter :: Rsearch = 2.,Rout=2.5 Real :: Xc,Yc,Xcv,Ycv,SLPmin,Distance,Vortmax,vartemp, Vortemp Integer :: Ifind,Istatus ! Date and time and data file name Type(datetime) :: s_date,c_date,e_date,temp_date integer, parameter :: fu=111, fv=112, fw=113, ft=114,fh=115,fslp=116 ! File ID for output integer :: irec1, irec2, orec,tstep !------------------------------------- ! Define data for holding the output Integer, parameter :: max_tc = 1000, max_obs=200, max_tc_ts=20 Type tcobs integer :: y,m,d,h,tc_num real :: lat, lon real :: vort, ocs, tano, TCoreano, slp, slpano End Type tcobs Type(tcobs), dimension(max_tc,max_obs) :: TCs ! This is the array for holding all the output! Type(tcobs) :: current_Obs, cur_Tc, last_TC Type(tcobs), dimension(max_tc_ts) :: current_tcs, predict_tcs ! 2 special array for holding TCs of current, prediction (+6h) integer, dimension(max_tc) :: num_obs ! Number of observation for each tc character(len=20),dimension(max_tc) :: tc_names integer :: num_tc ! Number of over all TC integer :: cur_num, pre_num, & ! Number of tc of current and previous time slice i,j,k, i_num,i_obs,i_tc ! Iterators logical :: iswarmcore, isstrong,isdvdz, DEBUG !!!!!!!!!!!!!!!!!!!!!!!!!!1 ! New version's variables ! real, parameter :: dx_search = 2.5, & ! Searching resoluton in Degree env_radius = 5. ! Half domain width to calc. enviromental fields integer :: di_search, & ! searching index resolution i1,i2,j1,j2, & ! Indices for calc. enviromental fields ic,jc, & ! Indices of analysing center env_irad, N_env,Nxx,Nyy ! real :: Slp_env,T_env logical :: exists !!!!! Kim Nguyen include 'netcdf.inc' integer ncid integer timest(2),timeco(2),timerco(2) real pi, rearth, r real, allocatable :: level(:),times(:) data pi/3.1415926536/,rearth/6371.22e3/,r/287./ integer lonid, latid,levid,timeid,nlon,nlat,nlevs,ntimes, & status,ierr,ier,mode,narch,iarch integer start(2),count(2) character*150 ifile data start/1,1/ !!!!!!!!!!!!!!!!!!!!!!!!! End !Namelist namelist /Datetimes/ s_date,e_date,ifile ! namelist /Domains/ Nx,Rx0,Dx,Ny,Ryfile,input_dir namelist /Criteria/ VorCrit,TCorecrit,H500crit,ocscrit,SLPcrit,CPradius,Cradius,DEBUG ! ! - Start program ! !--Read namelist information open(1,file="info.tcdetect") read(1,Datetimes) ! read(1,Domains) read(1,Criteria) close(1) ! ! open input netCDF outfile ! write(*,*)'file= ',ifile status=nf_open(ifile, nf_nowrite, ncid) if(status.ne.0) then write (*,*)' cannot open netCDF file; error code ',status write (*,*)' cannot open netCDF file; ',ifile stop end if write(*,*) "11111111TCdetect: Allocate fields..." ! turn OFF fatal netcdf errors call ncpopt(0) ! -------------------------------------- ! Get dimensions (co-ordinate variables) ! -------------------------------------- ! Pressure levels status = nf_inq_dimid(ncid, 'plev', levid) status = nf_inq_dimlen(ncid, levid, nlevs) status = nf_inq_varid(ncid, 'plev', levid) allocate (level(nlevs)) status = nf_get_vara_real(ncid, levid, 1, nlevs, level) write(*,*) level ! Longitudes status = nf_inq_dimid(ncid, 'lon', lonid) status = nf_inq_dimlen(ncid, lonid, Nx) status = nf_inq_varid(ncid, 'lon', lonid) allocate (Rx(Nx)) status = nf_get_vara_real(ncid, lonid, 1, Nx, Rx) write(*,*) Nx ! Latitudes status = nf_inq_dimid(ncid, 'lat', latid) status = nf_inq_dimlen(ncid, latid, Ny) status = nf_inq_varid(ncid, 'lat', latid) allocate (Ry(Ny)) status = nf_get_vara_real(ncid, latid, 1, Ny, Ry) ! Date and Time status = nf_inq_dimid(ncid, 'time', timeid) status = nf_inq_dimlen(ncid, timeid, ntimes) allocate (times(ntimes)) status = nf_get_vara_real(ncid, latid, 1, ntimes, times) narch = ntimes !write (*,*)Rx !write (*,*)Ry !:wallocate(Rx(Nx)); allocate(Ry(Ny)) allocate(Up(Nx,Ny,Np));allocate(Vp(Nx,Ny,Np));allocate(Wp(Nx,Ny,Np));allocate(Tp(Nx,Ny,Np)); allocate(Hp(Nx,Ny,Np));allocate(Vorp(Nx,Ny,Np));allocate(pTp(Nx,Ny,Np)); allocate(Slp1(Nx,Ny));allocate(Slp2(Nx,Ny));allocate(Vtang(Nx,Ny)); ! !-- Loop through times and read data on PRESSURE levels !-- and detect TCs ! c_date = s_date ! Current date tstep = 1 ! Time step num_tc = 0 ! Number of TCs pre_num = 0 ! Previous Number of TCs num_obs=0 ! Number of Observations !Initialize some pre-processing Dx = Rx(2)-Rx(1) di_search = dx_search/Dx env_irad = env_radius/Dx write(*,*) di_search write(*,*)"TCdetect: Start detecting..." iarch = 1 do while(date_comp(c_date,e_date) /=1 ) call histrd4(ncid,iarch,Nx,Ny,Np,'t',Nx,Ny,ptp) call histrd4(ncid,iarch,Nx,Ny,Np,'u',Nx,Ny,up) call histrd4(ncid,iarch,Nx,Ny,Np,'v',Nx,Ny,vp) call histrd4(ncid,iarch,Nx,Ny,Np,'hgt',Nx,Ny,hp) call histrd1(ncid,iarch,Nx,Ny,'mslp',Nx,Ny,slp1) iarch = iarch + 1 ! Calc relative vorticity do k=1,Np Call vorticity(Vorp(:,:,k),Up(:,:,k),Vp(:,:,k),Nx,Ny,Rx,Ry) Call Smooth_4p(Vorp(:,:,k),Nx,Ny) enddo !-- finish reading data !--- !-- Now, Time for TC detection! !--- write(*,*) write(*,'(" ------Timeslice(",i4,"), Date:",i4,"/",i2,"/",i2,":",i2)')tstep,c_date%y,c_date%m,c_date%d,c_date%h cur_num=0 loopi: do i=di_search,Nx-di_search,di_search loopj: do j=di_search,Ny-di_search,di_search if (DEBUG) write(*,*)"Searching at: ",i,j write(*,'(A)',advance="no")"." ! Step1: Search for SLP centers ! using birational interpolation, the SLP center is not neccessary one of the grid points. ! Radius of searching is 250km, If a SLP center is found, it serves as the center of TC !-Calculate domain indices for calculate environental domain i1 = i - env_irad i2 = i + env_irad j1 = j - env_irad j2 = j + env_irad if (i1<1) i1 = 1 if (i2>Nx) i2 = Nx if (j1<1) j1 = 1 if (j2>Ny) j2 = Ny Nxx=i2-i1+1 Nyy=j2-j1+1 ic=(i1+i2)/2 jc=(j1+j2)/2 call Search_center(-Slp1(i1:i2,j1:j2),Nxx,Nyy,Rx(i1:i2),Ry(j1:j2),Rx(ic),Ry(jc),Xc,Yc,SLPmin,Ifind,Rsearch) SLPmin=-Slpmin ! A slp center must be found if (0.eq.Ifind) then if (DEBUG) write(*,*)"!!! Cannot find SLP center" cycle loopj endif ! Determin the nearest grid point call X_to_i(Xc,ic,Rx,Nx) call X_to_i(Yc,jc,Ry,Ny) if (DEBUG) write(*,'(" SLP center: Lon:",f6.2,", Lat:",f6.2)') Rx(i),Ry(j) if (DEBUG) write(*,'(" SLP center: Lon:",f6.2,", Lat:",f6.2)') Xc,Yc if (DEBUG) write(*,'(" SLP center: i:",i6,", j:",i6)') ic,jc !-Calculate domain indices for calculate environental domain i1 = ic - env_irad i2 = ic + env_irad j1 = jc - env_irad j2 = jc + env_irad if (i1<1) i1 = 1 if (i2>Nx) i2 = Nx if (j1<1) j1 = 1 if (j2>Ny) j2 = Ny N_env=(i2-i1+1)*(j2-j1+1) if (DEBUG) write(*,*)i1,i2,j1,j2,N_env !-- SLP anomally must be at least 1hPa Slp_env = sum(Slp1(i1:i2,j1:j2)/N_env) SlpAno = SLPmin - Slp_env if (DEBUG) write(*,*)' SLP-anomaly',SLPmin,Slp_env,SlpAno ! Check if the center is already analysed? ! write(*,*)"Center:",Xc,Yc do i_tc = 1, cur_num cur_tc = current_tcs(i_tc) call geodistance(Xc,Yc,cur_tc%lon,cur_tc%lat,Distance) if (Distance .lt. Cradius) then if (DEBUG) then write(*,*)Xc,Yc write(*,*)"->>",Distance write(*,*) "!!! Already analysed for TC#",cur_tc%tc_num endif cycle loopj endif enddo if (SLPano>SLPCrit)then if (DEBUG) write(*,*)"!!! SLP-anomaly is too small" cycle loopj endif ! C-4: Check the variation of Tangt. wind according to height !--Caculate tangential wind around the center if (DEBUG) write(*,*)" Outer core streng (OCS) at 850mb" OCS=0 call calc_OCS(Up(i1:i2,j1:j2,1),Vp(i1:i2,j1:j2,1),Nxx,Nyy,Rx(i1:i2),Ry(j1:j2),Xc,Yc,OCS) if (DEBUG) write(*,*)Plev(k),OCS if (OCS < OCScrit) then if (DEBUG) write(*,*) "OCS is not strong enough",OCS cycle loopj endif ! Calculate Temperature pertubations and Tangential windstrength !Tano=0. !TCoreano= 1. ! Temporaty do k=1,Np Tanop(k) = ptp(ic,jc,k) - sum(ptp(i1:i2,j1:j2,k)/N_env) enddo TCoreano = .1*Tanop(1) + .2*Tanop(2)+.3*Tanop(3)+.4*Tanop(4) if (DEBUG) then write(*,*)' TCoreano',TCoreano endif ! C-3: Check warm core structure if (TCoreano < TCorecrit) then if (DEBUG) write(*,*)"!!! Core Temperature anomally is too small",TCoreano,TCorecrit cycle loopj endif ! TC-4 Now, a TC is indentified, is it a newly formed or an old one? current_Obs%lon=Xc current_Obs%lat=Yc current_Obs%slpano=slpano current_Obs%OCS=OCS current_Obs%tano=tano current_Obs%TCoreano=TCoreano current_Obs%slp=slpmin current_Obs%vort=Vortmax current_Obs%tc_num=0 ! Unknown ID yet! ! Check if the detected TC is a new one or a pre-existed one do i_tc = 1,pre_num call geodistance(Xc,Yc,predict_tcs(i_tc)%lon,predict_tcs(i_tc)%lat,Distance) if (Distance < CPRadius) then current_Obs%tc_num = predict_tcs(i_tc)%tc_num !--assign the ID of an old TC endif enddo !--new born TC if (current_Obs%tc_num.eq.0)then cur_num=cur_num+1 num_tc=num_tc+1 current_Obs%tc_num = num_tc current_tcs(cur_num)=current_Obs ! add new TC else !Check if is there a record of the same ID exists=.false. do i_tc = 1,cur_num if (current_Obs%tc_num .eq. current_tcs(i_tc)%tc_num) then exists=.true. ! choose the lower SPL TC and exit (without increasing current number of TC_obs) if (current_Obs%slp .lt.current_tcs(i_tc)%slp) then current_tcs(i_tc) = current_Obs endif exit endif enddo if (.not. exists) then cur_num=cur_num+1 ! increase current number of TC_obs -- only if not close the any TC! current_tcs(cur_num) = current_Obs ! add new TC endif endif write(*,*) write(6,'(" ----> Found TC #",i3," [",f5.1,":",f5.1,"]",f6.1)'),current_tcs(cur_num)%tc_num,Xc,Yc,SlpAno enddo loopj enddo loopi ! Update the global TCs array do i_tc = 1,cur_num cur_tc = current_tcs(i_tc) cur_tc%d=c_date%d cur_tc%m=c_date%m cur_tc%y=c_date%y cur_tc%h=c_date%h i = cur_tc%tc_num num_obs(i)=num_obs(i)+1 TCs(i,num_obs(i)) = cur_tc enddo ! Predict future tcs (update future_tcs) predict_tcs = current_tcs do i_tc = 1,cur_num cur_tc = current_tcs(i_tc) i = cur_tc%tc_num if (num_obs(i)>1) then ! If Number of observation is >=2, we predict the furture position ! By simple linear trend last_TC = TCs(i,num_obs(i)-1) predict_tcs(i_tc)%lon = 2*current_tcs(i_tc)%lon - last_TC%lon predict_tcs(i_tc)%lat = 2*current_tcs(i_tc)%lat - last_TC%lat endif enddo pre_num = cur_num cur_num = 0 ! if ((pre_num>0).and.(tstep>20)) then ! pause ! endif !--- !-- Finished TC detection! !--- call add_hrs(c_date,6) tstep=tstep+1 enddo status=nf_close(ncid) ! close(11) ! Finally, write to output file open(1, file="TCs.txt",status="unknown") write(1,"('TC#;',i5)")num_tc do i_num=1,num_tc write(1,*)"----------------------------------" write(1,"('TC#;',i5)")i_num write(1,"('Num_Obs:',i5)")num_obs(i_num) do i_obs=1,num_obs(i_num) cur_tc=TCs(i_num,i_obs) write(1,'(i4,";",i2,";",i2,";",i2,5(";",f9.2))')cur_tc%y,cur_tc%m,cur_tc%d,cur_tc%h,cur_tc%lon,cur_tc%lat, & cur_tc%TCoreano,cur_tc%slpano,cur_tc%ocs enddo enddo close(1) write(*,*) "TCdetect: Finished!" Stop end !********************************************************************* !read Netcdf data subroutine histrd1(histid,iarch,il,jl,name,ix,iy,var) character name*(*) integer start(3),count(3) real var(il,jl) ix=il iy=jl start(1) = 1 start(2) = 1 start(3) = iarch count(1) = ix count(2) = iy count(3) = 1 ! read data id = ncvid(histid,name,ierr) call ncvgt(histid,id,start,count,var,ierr) return ! histrd1 end !*************************************************************************** subroutine histrd4(histid,iarch,il,jl,kl,name,ix,iy,var) character name*(*) integer start(4),count(4),n real subvar(il,jl) real var(il,jl,kl) start(1) = 1 start(2) = 1 start(3) = 1 start(4) = iarch count(1) = il count(2) = jl count(3) = 1 count(4) = 1 ! read data id = ncvid(histid,name,ierr) do n=1,kl start(3)=n k=n call ncvgt(histid,id,start,count,subvar,ierr) do j=1,jl do i=1,il var(i,j,k)=subvar(i,j) end do end do end do return ! histrd4 end !********************************************************************* ! compute vorticity subroutine vort5(nx,ny,u,v,dx,dy,vort) real u(nx,ny),v(nx,ny),dx(nx,ny),dy(nx,ny),vort(nx,ny) ! ! calculate vorticities according to fourth-order accurate method ! from P.582 of "Mesoscale Meteorology and Forecasting", ed Ray. ! 222 format(a10,10(10(1x,f6.2),/)) do j=3,ny-2 do i=3,nx-2 deltax = dx(i,j) deltay = dy(i,j) dudy1 = 2.*(u(i,j+1) - u(i,j-1))/(3.*deltay) dudy2 = (u(i,j+2) - u(i,j-2))/(12.*deltay) dudy = dudy1 - dudy2 dvdx1 = 2.*(v(i+1,j) - v(i-1,j))/(3.*deltax) dvdx2 = (v(i+2,j) - v(i-2,j))/(12.*deltax) dvdx = dvdx1 - dvdx2 vort(i,j) = dvdx - dudy end do end do ! print 223,'vort5:vort= ',((vort(i,j),i=30,40),j=25,35) 223 format(a10,15(6(1x,E12.5),/)) return end ! Subroutine xy_to_ij using binary search to determine the gridpoint i nearest to X subroutine X_to_i(X,i,Rx,Nx) integer :: i,i1,i2,Nx real :: X real,dimension(Nx) :: Rx !Determin i i1=1 i2=Nx do while(i2-i1>1) i=(i1+i2)/2 if (X<=Rx(i)) then i2=i else i1=i endif enddo if (X-Rx(i1) < Rx(i2)-X ) then i=i1 else i=i2 endif end subroutine X_to_i