! ! Topical Cyclone Detection Program for regional climate models ! (C) Bui Hoang Hai, 2010 use datedef Implicit none !Dimension defines Integer :: Nx, Ny, Np Integer :: k850,k700,k500,k300 ! Integer, parameter :: Np=4 ! Real :: plev(Np) ! Data plev /850.,700.,500.,300./ Real, allocatable, dimension(:) :: Rx, Ry, Plev, Tanop, OCSp 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(:,:) :: slp,Slp2, Vtang,Vt,us,vs ! Variables for criteria and threholds !Real, dimension(Np) :: Tanop, OCSp Real :: Tcenter,Hcenter,Pcenter,Tout,Pout,Hout,SLPano,Tano,TCoreano,h500ano,OCS,Vtmax,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 VtCrit, & ! 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 = 5,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 = 3000, max_obs=2000, max_tc_ts=2000 Type tcobs integer :: y,m,d,h,tc_num real :: lat, lon real :: vort, ocs, tano, TCoreano, slp, slpano,Vtmax End Type tcobs Type(tcobs), dimension(max_tc,max_obs) :: TCs ! This is the array for holding all the output! Type(tcobs) :: 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 !!!!! 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 Kim Nguyen !Namelist namelist /Datetimes/ s_date,e_date,ifile namelist /Domains/ Nx,Rx0,Dx,Ny,Ryfile,input_dir namelist /Criteria/ VorCrit,TCorecrit,H500crit,VtCrit,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, 'lev', levid) status = nf_inq_dimlen(ncid, levid, nlevs) status = nf_inq_varid(ncid, 'lev', levid) allocate (Plev(nlevs)) allocate (Tanop(nlevs)) allocate (OCSp(nlevs)) status = nf_get_vara_real(ncid, levid, 1, nlevs, Plev) write(*,*) Plev Np = nlevs do k=1,Np if (Plev(k) .eq. 850. ) k850=k if (Plev(k) .eq. 700. ) k700=k if (Plev(k) .eq. 500. ) k500=k if (Plev(k) .eq. 300. ) k300=k enddo ! 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(slp(Nx,Ny));allocate(Slp2(Nx,Ny));allocate(Vtang(Nx,Ny));allocate(Vt(Nx,Ny));allocate(Us(Nx,Ny));allocate(Vs(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 write(*,*)"TCdetect: Start detecting..." iarch = 1 do while(date_comp(c_date,e_date) /=1 ) call histrd4(ncid,iarch,Nx,Ny,Np,'temp',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,'zg',Nx,Ny,hp) call histrd1(ncid,iarch,Nx,Ny,'psl',Nx,Ny,slp) call histrd1(ncid,iarch,Nx,Ny,'us',Nx,Ny,us) call histrd1(ncid,iarch,Nx,Ny,'vs',Nx,Ny,vs) 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 Vt (:,:) = sqrt(Us(:,:)**2 + Vs(:,:)**2) !-- finish reading data !--- !-- Now, Time for TC detection! !--- 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=2,Nx-1 loopj: do j=2,Ny-1 ! C-1: First Check first the maximum vorticity (850hPa) for the TC candidate ! Because of Hi-res, the criteria of Vorticity is higher (1.0e-4 instead of 1.0e-5) Vortemp = vorp(i,j,k850) if ( (vortemp > VorCrit) .and. & (vortemp >= vorp(i-1,j,k850)) .and. (Vortemp >= vorp(i+1,j,k850)) .and. & (vortemp >= vorp(i,j-1,k850)) .and. (Vortemp >= vorp(i,j+1,k850)) .and. & (vortemp >= vorp(i-1,j-1,k850)) .and. (Vortemp >= vorp(i+1,j+1,k850)) .and. & (vortemp >= vorp(i-1,j+1,k850)) .and. (Vortemp >= vorp(i+1,j-1,k850)) & ) then if (DEBUG) then write(*,*)"Candidate:",i,j,vorp(i,j,k850) write(*,'(" FG center :",f6.2,", Lat:",f6.2)') Rx(i),Ry(j) endif ! C-2: From the candidate point as first gess center, search for SLP center ! 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 call Search_center(-slp,Nx,Ny,Rx,Ry,Rx(i),Ry(j),Xc,Yc,SLPmin,Ifind,Rsearch) ! A slp center must be found if (0.eq.Ifind) then if (DEBUG) write(*,*)"!!! Cannot find SLP center" cycle loopj endif if (DEBUG) write(*,'(" SLP center: Lon:",f6.2,", Lat:",f6.2)') Xc,Yc call geodistance(Xc,Yc,Rx(i),Ry(j),Distance) if (Distance > 150.) then if (DEBUG) write(*,*)" SLP center and FG Center are to far away!" cycle loopj endif ! Check if the candidate is too close to boundary if (Xc+2.5>Rx(Nx) .or. Xc-2.5Ry(Ny) .or. Yc-2.5SLPCrit)then if (DEBUG) write(*,*)"!!! SLP-anomaly is too small" cycle loopj endif ! Calculate Temperature pertubations and Tangential windstrength Tano=0. if (DEBUG) write(*,*)' T-anomaly ' do k=1,Np call azimean(ptp(:,:,k),Nx,Ny,Rx,Ry,Xc,Yc,Rout,Tout,Tcenter,Tanop(k),iStatus) if (iStatus.eq.0) cycle loopj if (DEBUG) write(*,*)plev(k),Tanop(k) enddo Tano=Tanop(k850)+Tanop(k700)+Tanop(k500)+Tanop(k300) Tano=Tano/4 TCoreano= 0.1*Tanop(k850) + 0.2*Tanop(k700) + 0.3*Tanop(k500) + 0.4*Tanop(k300) if (DEBUG) then write(*,*)' T-anomaly',Tano write(*,*)' TCoreano',TCoreano endif ! calculate H anomaly at 500hpa call azimean(Hp(:,:,k500),Nx,Ny,Rx,Ry,Xc,Yc,Rout,Hout,Hcenter,H500ano,iStatus) if (iStatus.eq.0) cycle loopj if (DEBUG) write(*,*)' H500ano',H500ano,H500crit if (H500ano > H500crit) then if (DEBUG) write(*,*)"!!! Height anomaly at middle level is too small",H500ano,H500crit cycle loopj 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 !Check maximum wind Vtmax ! call Search_center(Vt,Nx,Ny,Rx,Ry,Xc,Yc,Xcv,Ycv,Vtmax,Ifind,Rsearch) ! if (DEBUG) write (*,*)'Max wind speed',Vtmax ! if (Vtmax < VtCrit) then ! if (DEBUG) write(*,*) " Wind strength is too small",Vtmax,VtCrit ! cycle loopj ! endif ! TC-4 Now, a TC is indentified, is it a newly formed or an old one? cur_num=cur_num+1 current_tcs(cur_num)%lon=Xc current_tcs(cur_num)%lat=Yc current_tcs(cur_num)%tc_num=0 current_tcs(cur_num)%slpano=slpano current_tcs(cur_num)%OCS=OCS current_tcs(cur_num)%tano=tano current_tcs(cur_num)%TCoreano=TCoreano current_tcs(cur_num)%slp=slpmin current_tcs(cur_num)%vort=Vortmax current_tcs(cur_num)%Vtmax=Vtmax ! 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 !--old TC current_tcs(cur_num)%tc_num = predict_tcs(i_tc)%tc_num endif enddo !--new born TC if (current_tcs(cur_num)%tc_num.eq.0)then num_tc=num_tc+1 current_tcs(cur_num)%tc_num = num_tc endif write(6,'(" ----> Found TC #",i3," [",f5.1,":",f5.1,"]")'),current_tcs(cur_num)%tc_num,Xc,Yc endif ! end of candidate loop 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(";",f22.2))')cur_tc%y,cur_tc%m,cur_tc%d,cur_tc%h,cur_tc%lon,cur_tc%lat, & cur_tc%tano,cur_tc%slpano,cur_tc%Vtmax 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